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Abstract 

To  support  the  Air  Force’s  Global  Reach  concept,  a  Common  Aero  Vehicle  is 
being  designed  to  support  the  Global  Strike  mission.  “Waypoints”  are  specihed  for 
reconnaissance  or  multiple  payload  deployments  and  “no-fly  zones”  are  specihed  for 
geopolitical  restrictions  or  threat  avoidance.  Due  to  time  critical  targets  and  mul¬ 
tiple  scenario  analysis,  an  autonomous  solution  is  preferred  over  a  time-intensive, 
manually  iterative  one.  Thus,  a  real-time  or  near  real-time  autonomous  trajectory 
optimization  technique  is  presented  to  minimize  the  flight  time,  satisfy  terminal  and 
intermediate  constraints,  and  remain  within  the  specihed  vehicle  heating  and  control 
limitations.  This  research  uses  the  Hypersonic  Cruise  Vehicle  (HCV)  as  a  simpli¬ 
fied  two-dimensional  platform  to  compare  multiple  solution  techniques.  The  solution 
techniques  include  a  unique  geometric  approach  developed  herein,  a  derived  analyti¬ 
cal  dynamic  optimization  technique,  and  a  rapidly  emerging  collocation  numerical  ap¬ 
proach.  This  up-and-coming  numerical  technique  is  a  direct  solution  method  involving 
discretization  then  dualization,  with  pseudospectral  methods  and  nonlinear  program¬ 
ming  used  to  converge  to  the  optimal  solution.  This  numerical  approach  is  applied  to 
the  Common  Aero  Vehicle  (CAV)  as  the  test  platform  for  the  full  three-dimensional 
reentry  trajectory  optimization  problem.  The  culmination  of  this  research  is  the  ver¬ 
ification  of  the  optimality  of  this  proposed  numerical  technique,  as  shown  for  both 
the  two-dimensional  and  three-dimensional  models.  Additionally,  user  implementa¬ 
tion  strategies  are  presented  to  improve  accuracy  and  enhance  solution  convergence. 
Thus,  the  contributions  of  this  research  are  the  geometric  approach,  the  user  imple¬ 
mentation  strategies,  and  the  determination  and  verification  of  a  numerical  solution 
technique  for  the  optimal  reentry  trajectory  problem  that  minimizes  time  to  target 
while  satisfying  vehicle  dynamics  and  control  limitation,  and  heating,  waypoint,  and 
no-fly  zone  constraints. 
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Common  Aero  Vehicle 


Autonomous  Reentry  Trajectory  Optimization 
Satisfying  Waypoint  and  No-Fly 
Zone  Constraints 

I.  Introduction 

1.1  Motivation 

Global  Strike  (GS)  and  Global  Persistent  Attack  (GPA)  are  two  of  the  seven 
United  States  Air  Force  Concepts  of  Operations  [104:pg.  10] .  Various  hypersonic  and 
reentry  vehicle  technologies  are  being  pursued  to  enable  such  a  prompt  global  reach  ca¬ 
pability  [118].  The  government  competition  for  the  Force  Application  Launch  from  the 
Continental  United  States  (FALCON)  [19:pg.4]  program  and  the  National  Aeronautics 
and  Space  Administration  (NASA)  Next  Generation  Launch  Technology  (NGLT)  Pro¬ 
gram  Office  [42]  validate  both  the  interest  and  need  for  future  research.  The  Common 
Aero  Vehicle  (CAV)  is  an  Unmanned  Aerial  Vehicle  (UAV)  supporting  these  technolo¬ 
gies  with  the  role  of  “Striking  from  Space”  [76].  One  goal  is  a  global  strike  mission 
where  targets  within  a  9,000  nautical  mile  range  can  be  reached  in  less  than  two 
hours  with  12,000  lb  payload  capabilities  [121].  To  analyze  integration  and  mission 
effectiveness,  the  Control  Sciences  Division  of  the  Air  Force  Research  Laboratory/Air 
Vehicles  Directorate  (AFRL/VA),  Aerospace  Vehicles  Technology  Assessment  and 
Simulation  Branch  at  Wright-Patterson  AFB  is  developing  the  Space  Access  Vehicles 
Mission  and  Operations  Simulation  (SAVMOS)  [47;  117].  SAVMOS  is  a  computer 
simulation  environment  designed  for  modeling  a  Military  Spaceplane  (MSP),  Com¬ 
mon  Aero  Vehicle  (CAV)  [57:pg.7;  32],  and  its  operations  system.  The  core  trajectory 
generation  for  reentry  is  the  Integrated  Development  and  Operations  System  (IDOS) 
[57;  47:pg.3].  This  allows  for  simulation  runs  from  reentry  to  the  target,  with  the 
trajectory  satisfying  vehicle  dynamics  and  heat /structural  load  constraints.  The  cur- 
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rent  CAV  implementation  within  SAVMOS  does  not  allow  for  rapid  mission/target 
modifications,  nor  is  this  a  current  operational  capability  for  Global  Strike. 

1.2  Problem  Description 

The  desire  to  compute  a  survivable  reentry  trajectory  has  been  around  since 
the  dawn  of  manned  space  flight.  The  Space  Shuttle  introduced  a  more  difficult 
challenge  by  including  a  larger  lift  capacity  reentry  vehicle  with  mission  needs  to 
precisely  navigate  to  a  desired  landing  site,  thus  requiring  appropriate  guidance  [39]. 
Since  then,  autonomous  reentry  vehicles  such  as  the  X-33  and  X-37  required  both 
autonomy  and  the  ability  to  accommodate  off-nominal  reentry  conditions  produced 
from  modeling  errors  or  vehicle  faults.  A  summary  of  the  on-going  reentry  trajectory 
methods  is  featured  in  [22],  One  article  of  interest  [36]  compared  five  of  the  pro¬ 
posed  trajectory  generation  techniques  and  tested  them  with  the  high-fidelity  X-33 
Marshall  Aerospace  Vehicles  Representation  in  C  (MAVERIC)  simulation  [38].  These 
results  were  re-presented  two  years  later  [37].  Clearly  the  civilian  objective  of  reen¬ 
try  through  non-hostile  environments,  with  the  single  mission  of  vehicle  recovery,  has 
been  addressed. 

For  military  applications  to  vehicles  such  as  CAV,  additional  requirements  are 
introduced.  First  of  all,  geopolitical  limitations  may  be  imposed,  thus  eliminating 
feasible  trajectories  that  would  overfly  these  restricted  areas.  Also,  reconnaissance  or 
multiple  payload  deployments  may  dictate  flying  through  specified  waypoints.  Way- 
points  may  also  be  used  to  adjust  the  reentry  trajectory  to  coordinate  with  other 
military  objectives  or  air  traffic  avoidance.  The  target  is  assumed  to  be  time-critical; 
therefore,  achieving  these  objectives  in  minimum  time  is  required. 

Therefore,  the  desired  capability  and  the  overall  goal  of  this  research  is: 

Given  a  time-sensitive  target,  impose  multiple  intermediate  waypoints,  add 
threat  avoidance  or  7io-fly  zones,  and  then  autonomously  generate  a  mini¬ 
mum  time  fly  able  and  survivable  trajectory  for  a  hypersonic  reentry  vehicle 
such  as  the  Common  Aero  Vehicle. 
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The  final  research  objective  is  to  prove  the  optimality  of  the  solution  to  the  above 
problem.  The  methodology  will  be  to  pose  it  as  a  minimization  problem  having  an 
objective  function  with  equality  and  inequality  constraints.  The  objective,  or  cost, 
function  includes  time  to  target;  however,  it  could  also  be  any  function  of  the  final 
states,  e.g.  terminal  velocity  or  energy.  The  waypoints  and  no-fly  zones  are  considered 
rigid  constraints,  meaning  waypoint  exact  passage  must  be  achieved,  and  no-fly  zones 
shall  not  be  violated.  A  survivable  trajectory  imposes  a  maximum  allowable  spike 
heating  constraint.  Lastly,  the  optimality  of  a  numerically  generated  solution  must 
be  verified.  Some  methods  enforce  optimality  in  the  solution  process,  others  will  need 
to  be  verified  after  the  solution  is  generated.  This  is  addressed  in  Chapters  III  and 
IV. 


1.2.1  Terminology.  The  overall  mission  objective  is  to  fly  from  an  initial 
point  to  a  final / terminal  point  or  target,  in  minimum  time.  The  generic  start  and 
finish  point  of  a  trajectory  are  called  endpoints.  Waypoints  are  specified  interme¬ 
diate  coordinates  to  fly  over  to  satisfy  payload  delivery  or  reconnaissance  mission 
requirements.  The  vehicle  must  fly  directly  over  each  waypoint,  also  called  waypoint 
passage ;  however,  time,  altitude,  bank  angle,  flight  path  angle,  velocity,  and  heading 
are  not  constrained.  Some  control  law  solutions  are  called  bang-bang  or  bang-level- 
bang,  meaning  the  control  takes  on  the  value  of  minimum,  zero,  or  maximum.  For  a 
bang-bang  controller  for  bank  angle,  a  turnpoint  is  defined  as  a  discrete  point  where 
there  is  a  change  in  the  current  constant  control;  therefore,  the  vehicle  rolls  into  a  turn 
or  completes  a  turn  at  a  turnpoint.  There  is  no  predetermined  limit  on  the  number 
of  discrete  changes  in  control,  i.e.  turnpoints.  If  waypoint  passage  occurs  at  some 
point  within  a  constant  control  turn,  the  waypoint  and  turnpoints  are  not  coincident. 
A  no-fly  zone  is  a  region  with  a  boundary  that  the  vehicle  may  contact  but  must 
not  violate.  The  entire  trajectory  can  be  broken  up  into  sequential  legs,  segments,  or 
phases ;  these  terms  may  be  used  interchangeably.  The  breaks  between  phases  may 
occur  at  waypoints,  turnpoints,  or  other  ending  criteria  such  as  no-fly  zone  contact. 
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The  problem  starts  at  the  initial  time  t0-  Since  there  can  be  multiple  waypoints  in 
a  mission,  each  is  numbered  i,  and  the  time  of  passage  is  denoted  tt.  Similarly,  each 
no-fly  zone  is  numbered  j:  however,  the  time  of  interest  for  each  no-fly  zone  is  the 
time  of  no-fly  zone  boundary  contact  tr  The  target  is  reached  at  the  final  or  terminal 
time  tf  and  is  called  target  intercept. 

1.2.2  Hypersonic  Cruise  Vehicle  (HCV)  Overview.  The  Hypersonic  Cruise 
Vehicle  (HCV)  is  a  proposed  hypersonic  vehicle  to  meet  specified  global  reach  ob¬ 
jectives.  The  HCV  objective  is  to  deliver  12,000  lb  of  payload  at  a  range  of  9,000 
nautical  miles  (nmi)  in  2  hours  [20;  105],  at  a  notional  altitude  of  100,000  feet.  The 
HCV  has  thrust  thus  validating  a  simplifying  constant  altitude  assumption  as  used 
herein.  The  HCV  provides  a  platform  with  a  large  turn  radius,  making  it  ideal  to 
validate  optimal  trajectories  in  the  two-dimensional  horizontal  plane.  The  HCV  is  of 
additional  interest  because  it  may  have  the  mission  of  deploying  a  CAV  [121].  More 
details  of  the  HCV  specifications  and  limitations  are  provided  in  Section  3.3. 

1.2.3  Common  Aero  Vehicle  (CAV)  Overview.  The  following  abstract  pro¬ 
vides  a  succinct  introduction  to  the  Common  Aero  Vehicle: 


The  Common  Aero  Vehicle  (CAV)  is  a  concept  which  describes  a  space 
reentry  aeroshell  launched  into  space  on  a  suitable  vehicle,  which  then 
survives  atmospheric  reentry,  reduces  its  speed  to  low  Mach  numbers  or 
even  sub-Mach,  and  dispenses  a  cargo,  payload  or  weapon  in  the  Earth’s 
atmosphere.  The  conceptual  CAV  might  be  propelled  into  space  by  any  of 
a  number  of  present  or  future  launch  platforms  and  dispense  a  wide  variety 
of  cargoes,  payloads  or  weapons.  Development  of  the  CAV  capability  will 
satisfy  future  requirements  enunciated  in  numerous  national  visions,  future 
studies,  and  military  plans  and  will  ultimately  be  necessary  to  fully  realize 
the  opportunities  inherent  in  operating  from,  through,  and  in  space  and 
give  true  meaning  to  the  phrases  global  reach,  global  power  projection, 
and  global  engagement.  [75] 

1.2.4  Launch  Scenario.  The  launch  profile  of  study  is  a  sub-orbital  reentry. 

Launch  may  be  achieved  via  an  expendable  launch  vehicle,  expendable  air  launch 
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rocket,  the  notional  reusable  Space  Operations  Vehicle  (SOV)  [75;  86:pg.4],  or  the 
Hypersonic  Cruise  Vehicle  (HCV)  [121].  A  benefit  of  sub-orbital  launch  is  that  the 
launch  vehicle  can  be  selected  uniquely  based  on  the  range  to  the  target  or  other  mis¬ 
sion  dependent  variables.  Additionally,  development  can  be  performed  incrementally, 
meaning  a  lesser  performance  launch  vehicle  can  support  the  CAV  now,  and  later 
advancements  in  launch  vehicle  technologies  can  be  directly  applied  to  future  CAV 
missions.  HCV  trajectory  optimization  is  performed  as  a  simplified,  constant  alti¬ 
tude,  hypersonic  launch  vehicle.  This  analysis  is  applicable  to  a  launch  vehicle/CAV 
combined  mission,  as  computed  in  [121],  or  as  an  initial  guess  to  a  CAV  only  reentry 
trajectory. 

1.2.5  Computational  Objectives.  The  objective  of  this  research  effort  is 
to  provide  a  solution  using  the  computational  capacity  that  could  be  aboard  an  au¬ 
tonomous  reentry  vehicle  such  as  the  CAV.  The  current  belief  is  this  would  be  compa¬ 
rable  to  the  computing  power  of  a  desktop  personal  computer.  The  uplink  capability 
of  the  vehicle  would  then  only  require  receiving  mission  specific  details,  and  not  an 
entire  mission  solution. 

1.2. 5.1  Real-time  updates.  Even  though  the  time  from  separation  to 
target  may  be  minimal,  enroute  mission  updates  may  still  be  a  possibility.  Addi¬ 
tional  work  in  this  held  is  being  performed  and  is  termed  Real-Time  Optimal  Control 
(RTOC)  [10;  79].  Therefore,  rapid  onboard  computing  is  also  justified  for  real-time 
updates  synergistic  with  the  vehicle  inner-loop  control  system.  The  inner-loop  system 
is  assumed  to  have  the  capability  to  adjust  for  minor  deviations  due  to  wind  while 
maintaining  the  computed  trajectory  profile.  The  research  herein  is  compatible  with, 
but  does  not  address,  RTOC  and  inner-loop  guidance. 

1.2. 5. 2  Autonomous  Pre-conflict  Analysis.  Autonomy  and  minimal 
computation  power  are  key  elements  to  enabling  hypothetical  pre-conflict  scenario 
analysis.  A  wargaming  capability  is  a  key  element  of  the  Course  of  Action  (COA) 
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Analysis  phase  and  COA  Comparison  phase  of  the  Joint  Air  Estimate  Process  ( JAEP) 
[103:pg.III-4]  and  [102:pg.  1] .  This  research  provides  an  enabling  technology  towards 
providing  this  wargaming  objective. 

1.3  Technology  Advancement 

In  order  to  advance  current  technology  and  contribute  uniquely  to  the  engi¬ 
neering  held,  a  problem  or  scenario  must  be  identified  that  has  not  been  previously 
achieved  and  has  a  target  audience  or  customer.  Section  1.1  introduced  reentry  tra¬ 
jectory  generation  work  being  done  but  also  pointed  out  that  the  work  is  catered  to 
commercial  applications  and  neglects  unique  military  concerns  such  as  threat  regions 
or  geopolitical  no-fly  zones.  With  the  military  utility  of  the  Global  Positioning  Sys¬ 
tem  (GPS),  many  UAVs  such  as  Global  Hawk  and  Predator  [112:pg.2],  achieve  route 
definition  through  waypoints.  Various  trajectory  optimization  strategies  navigating 
through  waypoints  have  been  studied  [3;  46;  106;  113;  119]  and  are  addressed  in  Sec¬ 
tion  2.2.  The  simulation  environment  SAVMOS,  mentioned  in  Section  1.1,  integrates 
one  of  the  leading  trajectory  generation  software  packages  [17;  29],  but  fails  to  pro¬ 
vide  autonomous  trajectory  generation,  waypoint  satisfaction,  and  threat  avoidance. 
Therefore,  the  proposed  enabling  technology  is  to  develop  an  autonomous  trajectory 
generation  algorithm  to  satisfy  vehicle  dynamics,  survivability,  waypoint,  and  threat 
or  no-fly  zone  constraints  as  stated  in  Section  1.2.  Autonomy  allows  for  onboard  real¬ 
time  computations  as  well  as  military  capability  analysis,  i.e.  wargaming  as  addressed 
in  Section  1.2. 5. 2.  The  case  study  vehicles  for  the  work  proposed  herein  are  the  HCV 
and  the  CAV  described  in  Sections  1.2,  1.2.3,  and  1.2.4  due  to  the  HCV’s  cruise  capa¬ 
bility  and  the  CAV’s  global  strike  capability;  commensurate  with  the  research  interest 
expressed  by  AFRL,  Air  Vehicles  Directorate  (AFRL/VA)  and  AFRL,  Space  Vehicles 
Directorate  (AFRL /VS). 
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1.4  Document  Outline 

The  objectives  of  this  research  are  to  determine  and  verify  a  viable  optimal 
solution  technique  that  computes  the  minimum  time  reentry  trajectory  which  strikes 
the  target  while  satisfying  the  vehicle  dynamics,  heating  limitation,  intermediate  way- 
points,  and  no-fly  zone  constraints.  A  unique  aspects  of  this  research  is  the  combina¬ 
tion  of  these  multiple  constraints  into  one  solution,  and  verifying  the  optimality  of  a 
numerically  generated  solution.  Before  approaching  the  problem  as  a  whole,  each  one 
of  the  independent  components  is  researched.  The  previous  research  on  any  single 
task  is  summarized  in  Chapter  II.  In  addition  to  the  research  summary  is  a  description 
of  how  well  that  solution  technique  may  contribute  to  the  combined  problem.  The 
end  of  the  chapter  includes  the  decision  to  pursue  the  dynamic  optimization  solution 
technique  and  the  justification  of  its  application  to  this  problem.  Chapter  III  contains 
a  generic  dynamic  optimization  problem  definition,  followed  by  vehicle  and  mission 
specific  constraints.  The  problem  complexity  begins  with  a  simpler  two-dimensional 
model  and  then  increases  to  a  more  complex  three-dimensional  model.  The  simplest 
problem  setup  is  the  HCV  at  constant  altitude  and  constant  speed.  Some  complex¬ 
ity  is  then  introduced  by  allowing  the  speed  of  the  HCV  to  decrease  throughout  the 
trajectory.  The  CAV  is  the  final  vehicle  considered  with  altitude  no  longer  confined; 
thus  with  a  change  in  altitude  comes  the  inclusion  of  the  heating  constraint.  Chapter 
IV  contains  the  solution  development,  both  analytical  and  numerical.  In  spite  of  the 
initial  focus  on  the  dynamic  optimization  technique,  an  advantageous  2-D  geometric 
solution  is  also  derived  and  presented.  Following  this  contribution  is  the  full  3-D 
numerical  solution.  Subsequently,  the  optimality  of  these  numerical  results  is  veri¬ 
fied  using  the  previously  derived  analytical  criteria.  An  elaborate  problem  setup  is 
required  to  prove  optimality  to  the  level  of  detail  described  above;  however,  once  the 
numerical  solution  process  is  verified,  a  simpler  setup  is  implemented.  Before  finishing 
this  chapter,  the  numerical  results  from  a  streamlined,  user  representative,  problem 
setup  is  compared  to  the  previous,  more  rigorous,  results.  Chapter  V  provides  the 
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conclusions  to  be  drawn  from  the  research  herein,  contributions,  and  comments  on 
future  research. 

The  appendices  provide  additional  details  to  complement  the  material  in  the 
main  body.  These  include  specifics  on  the  CAV  aerodynamic  model  (Appendix  A), 
nondimensionalization  of  the  equations  of  motion  (Appendix  B),  derivation  of  the 
shortest  path  between  points  (Appendix  C),  determination  of  the  geometrically  op¬ 
timal  waypoint  passage  point  (Appendix  D),  important  solution  aspects  of  the  pseu- 
dospectral  methods  (Appendix  E),  and  strategies  for  efficient  user  implementations 
of  the  numerical  solvers  (Appendix  F). 


II.  Previous  Research 


This  research  provides  a  synergistic  solution  encompassing  three  key  technolo¬ 
gies;  trajectory  generation,  waypoint  satisfaction,  and  threat  or  no-fly  zone  avoidance. 
The  following  addresses  the  research  in  each  area  separately  and  discusses  the  appli¬ 
cability  to  a  combined  solution. 

2.1  Reentry  Trajectory  Generation 

2.1.1  Space  Shuttle.  The  foundation  of  much  of  the  current  lifting  body 
reentry  vehicle  research  [24;  55;  61;  81;  84;  89;  90;  91;  92;  122]  is  based  on  or  at  least 
baselined  from  the  Space  Shuttle  [39].  The  Shuttle  trajectory  solution  was  based  on 
a  drag  acceleration  envelope  which  represented  the  various  vehicle  limitations;  such 
as  surface  temperature,  dynamic  pressure,  load  limits,  and  glide  limitations  as  seen 
in  Figure  1(a).  The  trajectory,  as  a  function  of  velocity,  is  broken  up  into  phases  such 
as  temperature  control,  equilibrium  glide,  and  constant  drag  to  traverse  through  the 
specified  envelope,  as  seen  in  Figure  1(b). 


RELATIVE  SPEED  THOUSANDS  OF  FPS 


(a)  Flight  Corridor 


(b)  Reference  Drag  Profile 


Figure  1:  Shuttle  Entry  -  Operational  Angle-of- Attack  Profile.  [39] 


2.1.2  X-33  and  X-37.  The  list  of  research  referenced  in  the  previous  para¬ 
graph  and  [82]  uses  the  X-33  as  the  test  vehicle  but  others  [1;  16;  28;  45;  60;  63] 
include  the  X-37  as  an  applicable  vehicle.  For  reusable  launch  vehicles  [83;  84], 
specifically  X-40A  [85],  fault  mitigation  is  another  pursued  topic.  One  paper  [37] 
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Figure  2:  X-33  -  Ground  Track  for  Six  Analyzed  Entries  [81] 


analyzes  multiple  reentry  trajectory  generation  techniques  and  provides  test  results 
comparing  five  methods  [24;  35;  81;  90;  122],  The  core  techniques  are  the  baseline 
guidance  method  [35]  like  the  Shuttle,  the  linear  quadratic  method  [24],  the  numerical 
predictor-corrector  method  [122],  the  drag-energy  three-dimensional  method  [81],  and 
the  quasi-equilibrium  glide  method  [90].  The  algorithms  producing  the  best  results 
were  the  quasi-equilibrium  glide  [90]  and  the  Evolved  Acceleration  Guidance  Logic 
for  Entry  (EAGLE)  [54;  81]  which  are  shown  in  Figure  2.  The  quasi-equilibrium 
glide,  having  apparent  roots  from  pseudo-equilibrium  glide  [39],  appears  to  be  very 
promising;  however,  the  EAGLE  program  has  had  more  extensive  testing  and  simu¬ 
lation  work  for  validation.  Furthermore,  this  is  the  algorithm  currently  incorporated 
in  SAVMOS  with  Mat  lab®  implementation.  The  Program  To  Optimize  Simulated 
Trajectories  II  (POSTII)  [71]  is  trajectory  optimization  and  simulation  tool  [21;  65]. 
It  has  a  two-point  boundary  value  problem  solver;  however,  there  is  not  a  readily 
available  means  to  incorporated  waypoints  or  no-fly  zone  constraints.  The  Optimal 
Trajectories  by  Implicit  Simulation  (OTIS)  [59;  67]  software  by  Glenn  Research  Center 
uses  the  newly  emerging  collocation  method  for  satisfying  the  differential  equations 
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Figure  3:  Fuzzy  Logic  -  Simulation  of  4  Waypoints  Trajectory  [3] 

and  solving  the  two-point  boundary  value  problem.  This  software  solves  for  a  set 
of  parameters  to  optimize  the  specified  trajectory  criteria.  One  critic  [62]  comments 
that  OTIS  does  not  enforce  continuity  over  phase  transitions;  however,  the  software 
package  used  herein  can  accept  user  defined  functions  to  remain  continuous  across  a 
phase  break.  This  is  described  more  in  Section  3.5  and  Appendix  F. 

2.1.3  Feedback  Linearization.  Feedback  linearization  transforms  a  nonlinear 
system  into  a  fully  or  partially  linear  system,  and  then  uses  linear  design  techniques 
to  complete  the  control  design  [93:pg.204,207] .  Feedback  linearization  can  also  be 
used  to  provide  an  entry  tracking  law  [9].  Feedback  linearization  is  most  applicable 
when  linearizing  about  a  nominal  point  or  trajectory;  however,  a  nominal  trajectory 
does  not  yet  exist.  Thus  this  approach  is  too  susceptible  to  the  non-linearities  of  this 
reentry  trajectory  problem. 

2. 2  Waypoints 

2.2.1  Fuzzy  Logic.  A  non-traditional  approach  to  waypoint  control  is  a 
five-dimensional  waypoint-based  fuzzy  guidance  system  (FGS)  for  unmanned  aircraft 
vehicles  [3;  46].  The  heading  and  altitude  for  each  waypoint  are  specified,  either  as  a 
restriction  or  an  enhancement.  Figure  3  shows  a  spiral  descent  was  correctly  identified 
when  two  consecutive  waypoints  were  otherwise  non-flyable.  This  approach  is  valid, 
but  the  other  facets  of  the  proposed  problem  can  not  readily  be  incorporated  into  an 
overarching  analytical  formulation. 
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2-D  Smooth  and  Optimal  Trajectory 
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Figure  4:  Riccati  Equation  -  Time  Optimal  Trajectory  of  the  Horizontal  Flight 

with  Constant  Speed  [119] 


2.2.2  Discontinuous  Lagrange  Multipliers.  A  discontinuous  Lagrange  multi¬ 
plier  is  one  that  has  a  jump,  or  discontinuity,  at  an  interior-point  constraint,  e.g. 
a  waypoint.  A  discontinuous  Lagrange  multiplier  arises  when  trying  to  solve  a 
two-point  boundary  problem  that  now  has  one  or  more  interior-point  constraints 
[14:pg.l01;  44:pg.275].  One  author  [119]  combined  discontinuous  Lagrange  multipli¬ 
ers  with  linearization  about  an  initial  guess  and  then  propagated  the  Riccati  differ¬ 
ential  equation  to  compute  components  for  the  Lagrange  multipliers.  This  approach 
is  appealing  to  the  problem  posed  herein  because  analytical  gradients  can  be  directly 
computed  since  a  method  for  computing  the  Lagrange  multipliers  exists.  This  ap¬ 
proach  is  used  for  the  two-dimensional  analytical  results  and  for  validation  of  the 
three-dimensional  numerical  results. 


2.2.3  Linearization  and  Propagating  Riccati  Equation.  Linearizing  the  sys¬ 
tem  allows  use  of  established  linear  optimization  and  control  techniques.  One  example 
is  establishing  line  segments  between  waypoints  and  then  designing  a  linear  quadratic 
regulator  (LQR)  for  line-following  guidance  [113].  Bryson  introduces  a  solution  by 
Riccati  equation  [12:pg.206]  and  expands  it  to  following  a  desired  output  [12:pg.217]. 
In  another  text  [13:pg.93],  Bryson  presents  the  algorithm  as  a  backward  information 
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filter-smoother  and  this  is  the  process  used  in  Figure  4  [119].  A  receding  horizon 
planning  strategy  using  Mixed-Integer  Linear  Programming  (MILP)  [106]  is  used  to 
command  a  UAV  through  waypoints  and  threat  zones.  This  appears  very  similar  to 
the  proposed  research;  however,  this  technique  takes  advantage  of  loiter  circles  and 
implements  intermediate  courses  based  on  limited  current  information.  The  linearity 
of  the  system  in  the  local  region  and  the  step  size  during  an  iteration  can  drastically 
affect  the  performance  of  a  linearization  technique;  therefore,  a  linearized  solution 
will  not  be  used  to  find  a  solution  for  this  research. 

2.2.4  Steepest  Descent.  Steepest  descent  is  one  of  the  simplest  gradient 
based  algorithms,  but  it  is  sensitive  to  scaling  parameters  [13: pg .28,125,1 54]  and  is 
notoriously  slow  to  converge  for  problems  with  long  narrow  valleys  [58:pg.3-5].  For 
these  reasons,  alternate  gradient  based  methods  are  pursued  for  this  research. 

2.2.5  Sequential  Quadratic  Programming.  Sequential  Quadratic  Program¬ 
ming  (SQP)  is  also  gradient  based,  and  SQP  methods  represent  the  state-of-the-art  in 
nonlinear  programming  methods  [58:pg.3-29].  By  computing  the  Hessian,  typically  in¬ 
directly  via  the  Broyden,  Fletcher,  Goldfarb,  and  Shanno  (BFGS)  method  [107:pg.293], 
second-order  gradient  information  is  used  that  dramatically  improves  performance 
over  the  first-order  steepest  decent  method.  As  an  example,  a  minimum-time-to- 
turn  flight  algorithm  is  computed  for  an  F-4  fighter  aircraft  using  SQP  [66].  Due  to 
the  fast  convergence,  this  is  the  core  algorithm  embedded  in  many  numerical  solu¬ 
tion  techniques.  For  example,  an  industry  standard  from  Stanford,  Sparse  Nonlinear 
Optimizer  (SNOPT),  is  based  on  SQP. 

2.3  No-Fly  Zones 

Threat  or  no-fly  zones  are  important  in  multiple  contexts.  The  military  threat 
zone  is  typically  understood  as  a  threat  severity  being  proportional  to  the  distance 
from  the  threat  [49;  56;  100;  108],  as  seen  in  Figure  5.  Other  instances  of  no-fly  zones 
or  regions  are  obstacle  avoidance  [120]  and  air  traffic  avoidance  [70].  Incorporation  of 
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Figure  5:  Threats  -  Ground  Track  for  Horizontal-Plane  Optimization  [108] 

the  no-fly  zones,  within  the  trajectory  planning  with  waypoints  is  the  subject  of  the 
research. 

2.3.1  Voronoi  Diagram.  A  threat  avoidance  grid  around  threat  sources  is 
first  segmented  into  equal  threat  distance  segments,  or  more  accurately  described  as 
a  Voronoi  Diagram  [49],  as  shown  in  Figure  6.  This  approach  allows  for  incrementally 
solving  the  solution,  meaning  a  local  solution  exists  before  completing  the  terminal 
phase  computation.  The  authors  [49]  use  a  spline  algorithm,  with  results  shown  in 
Figure  6(b),  to  convert  the  sharp  corners  of  the  Voronoi  diagram  solution  in  Figure 
6(a)  to  smooth  flyable  trajectories.  Algorithms  for  creating  these  Voronoi  diagrams 
are  available  in  Matlab®  [25;  64:pg.65].  This  works  well  for  threats,  modeled  as 
sources,  but  does  not  readily  incorporate  waypoints,  which  may  be  sinks  in  cost. 
Therefore,  this  technique  is  not  incorporated  within  this  research. 

2.3.2  Barrier  and  Interior  Point  Methods.  This  approach  may  be  used  for 
either  threat /terrain  avoidance  [100]  as  shown  in  Figure  7,  or  air  traffic  avoidance  [70] 
as  shown  in  Figure  8.  The  barrier  and  interior  point  methods  [5:pg.380]  increase 
the  cost  as  the  constraint  boundary  is  approached.  The  solution  process  starts  on 
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Figure  6:  Voronoi  Diagram  [49] 


Figure  7:  Barrier  -  Terrain  with  Threats  Formulated  as  an  Exponential  Function. 
On  the  Right  is  an  Overhead  View  [100]. 
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Figure  8:  Interior  Point  -  Projection  onto  yz  and  xy  Plane  of  the  Aircraft  Trajec¬ 
tories  when  No-Fly  Zone  is  Modeled  as  an  Ellipse  [70] 


the  interior,  i.e.  within  the  feasible  region,  and  attempts  to  improve  the  cost.  An 
exponential  or  1/x  function  is  typically  applied  that  approaches  infinity  at  the  con¬ 
straint  boundary;  therefore,  such  high  cost  at  the  boundary  drives  the  solution  to 
remain  within  the  interior.  A  downside  to  this  method  is  the  cost  function  is  either 
undefined  or  infinity  for  constraint  violations;  consequently,  incidental  iteration  steps 
within  the  constraint  could  cause  a  solution  algorithm  to  diverge.  A  similar  approach 
uses  a  penalty  function  that  increases  cost  for  violating  the  constraint,  for  example 
/  =  Aexp(— r2/b)  is  the  penalty  shown  in  Figure  7.  Here,  A  is  amplitude,  r  is  ra¬ 
dius,  and  b  is  a  scaling  factor  to  adjust  the  width.  Other  functions  such  as  a  step 
or  polynomials  [64:pg.441]  can  also  be  used  to  penalize  departure  from  the  feasible 
region,  i.e.  violating  the  constraint.  A  function  representing  the  constraint  that  has 
continuous  first  derivatives  [94:pg.99]  is  conducive  to  analytical  gradients.  The  an¬ 
alytical  development  herein  does  not  use  penalty  functions;  however,  such  weights 
for  feasible  violations  are  internal  to  the  numerical  algorithms  in  the  form  or  merit 
functions  [6;  8]. 


2.3.3  Receding  Horizon  and  Artificial  Intelligence.  Several  authors  [52;  53; 
106]  approach  the  UAV  using  obstacle  problems  with  a  receding  horizon  controller 
(RHC).  This  technique  is  based  on  the  most  current  information,  i.e.  waypoints  and 
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Figure  9:  Receding  Horizon  -  Representation  of  a  Cost-to-Go  Function  [53;  74] 
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Figure  10:  A*  Search-Optimal  Flight  Trajectory  in  Stationary  Obstacles  [120] 
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constraints,  which  may  only  extend  out  radially  to  a  limited  “detection  horizon”,  to 
find  a  solution  out  to  the  “planning  horizon”  shown  in  Figure  9.  The  use  of  an  admis¬ 
sible  set  of  allowed  solutions  allows  the  vehicle  to  negotiate/reject  feasible  solutions 
but  may  be  trapped  by  obstacles.  Ensuring  goal  attainment,  even  in  a  disturbance 
environment,  is  the  objective  of  these  planners.  The  shortest  path  from  each  node 
and  cost-to-go  are  similar  to  a  weighted  Voronoi  diagram.  This  avoidance  of  block  ob¬ 
stacles  seems  most  applicable  to  small  maneuverable  UAVs,  such  as  helicopters,  in  an 
obstacle  rich  urban  environment.  This  is  exactly  the  case  in  Figure  10  for  trajectory 
planning  using  an  A*  search  scheme  in  artificial  intelligence  [120].  For  the  research 
herein,  the  waypoints  and  threats  are  known  throughout  the  trajectory;  therefore, 
the  receding  horizon  approach  will  not  be  used  since  the  entire  trajectory  will  be 
computed  initially. 

2-4  Pseudo  spectral  Methods 

The  pseudospectral  methods  are  generic  numerical  solution  techniques,  not  spe¬ 
cific  to  any  one  class  of  problems  addressed  in  the  previous  sections,  ft’s  one  piece 
of  solving  the  optimal  control  problem  via  the  direct  method  instead  of  the  indirect 
method.  Two  authors  [12;  14;  44]  demonstrate  many  examples  of  the  indirect  ap¬ 
proach;  meaning  solving  a  series  of  optimality  conditions  to  then  indirectly  arrive  at 
the  optimal  control.  Except  for  very  simple  problems  that  have  a  purely  analytical 
solution,  discretization  will  certainly  be  part  of  the  solution  process,  e.g.  using  nu¬ 
merical  integration.  Defining  the  optimality  criteria  has  been  called  “dualization”  [27] 
since  it  creates  a  set  of  costates  or  dual  variables.  The  problem  arise  where  the  states 
are  known  at  the  initial  time,  and  the  costates  are  specified  at  the  final  time.  The 
shooting  method  is  one  indirect  approach  that  integrates  the  states  forward  and  the 
costates  backwards;  continuing  until  converging  to  a  solution  that  satisfies  the  initial 
and  final  conditions  for  the  states  and  costates,  respectively.  This  entails  Runge- 
Kutta  numerical  integration  which  is  a  form  of  discretization.  There  are  other  ways 
to  solve  the  problem  by  converting  the  unknown  quantities  of  the  indirect  optimiza- 
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tion  problem  into  a  parameter  optimization  problem  [43];  thus  solving  for  a  discrete 
number  of  variables.  In  summary,  the  indirect  method  specifies  a  set  of  optimality 
conditions  first,  then  discretizes  the  problem  to  find  the  optimal  control.  Conversely, 
the  direct  method  interchanges  the  dualization  and  discretization  steps  as  described 
next. 

There  are  at  least  five  authors  that  use  the  direct  approach  to  solve  the  op¬ 
timal  control  problem,  each  have  developed  their  own  software  package.  Professor 
O.  von  Stryk’s  papers  [95;  96]  discuss  the  direct  versus  indirect  method  and  provide 
example  problems  solved  with  his  direct  collocation  software  (DIRCOL).  DIRCOL 
is  also  used  in  [50]  to  solve  a  Space  Shuttle  heating  problem.  John  Betts  [6;  7;  8] 
also  converts  the  problem  into  a  nonlinear  program  (NLP)  and  created  the  Sparse 
Optimal  Control  Software  (SOCS)  as  a  solver.  The  software  package  called  Nonlinear 
Trajectory  Generation  (NTG)  [68]  focuses  on  solving  a  lower  dimensional  problem, 
also  satisfying  conditions  at  “knots”  or  collocation  points.  All  of  these  software  pack¬ 
ages  rely  on  the  NLP  solvers  Nonlinear  Program  Solver  (NPSOL)  or  SNOPT  from  the 
Stanford  Systems  Optimization  Laboratory  (SOL).  I.  Michael  Ross  [27;  77;  79]  and 
David  Benson  et  al.  [4]  use  the  Legendre-Gauss-Lobatto  (LGL)  and  Legendre- Gauss 
(LG)  pseudospectral  methods,  respectively.  The  pseudospectral  method  computes  the 
states,  controls,  and  derivatives  at  discrete  nodes;  then  uses  the  Karush-Kuhn- Tucker 
(KKT)  conditions  to  enforce  optimality.  The  benefit  to  the  pseudospectral  method 
is  that  the  solution  to  the  differential  equations  is  computed  simultaneously  across 
all  nodes  rather  than  sequentially  via  Runge-Kutta  integration;  therefore  alleviating 
the  potentially  debilitating  instability  issue  of  the  forward  and  backward  integration 
in  the  shooting  method.  These  authors’  respective  software  packages  are  DIDO  and 
Gauss  Pseudospectral  Optimal  Control  Software  (GPOCS)  [72],  One  of  the  challenges 
to  the  direct  method  is  obtaining  the  costate  information.  These  software  packages 
also  have  costate  information;  specifically,  DIDO  uses  the  Covector  Mapping  Theorem 
and  GPOCS  uses  Costate  Mapping  Theorem  to  convert  the  discretized  KKT  multi¬ 
pliers  to  the  continuous  Lagrange  multipliers,  see  Figure  11.  Figure  11  also  shows  the 
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similarity  between  these  latter  two  approaches.  The  Lagrange  multipliers  are  the  key 
to  verifying  optimality  used  within  this  research. 
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Figure  11:  Computation  of  the  Continuous  Lagrange  Multipliers 


2.5  Concluding  Remarks 

Various  methods  are  shown  to  solve  any  single  research  objective;  however,  the 
challenge  is  finding  a  method  to  combine  the  vehicle  dynamics,  heating,  waypoints, 
and  no-fly  zone  constraints.  Dynamic  optimization  can  combine  these  constraints  into 
a  single  problem,  but  many  numerical  solution  techniques  have  pitfalls  that  hinder 
convergence.  The  direct  method,  using  pseudospectral  methods,  avoids  some  of  the 
numerical  pitfalls,  and  thus  is  extremely  appealing  to  solving  the  optimal  control 
problem.  Obtaining  the  costate  information  from  the  proposed  direct  method  is 
critical  to  verifying  optimality.  Programs  such  as  DIRCOL,  DIDO,  and  GPOCS 
provide  such  costate  information.  The  challenge  within  this  research  is  to  formulate 
the  dynamic  optimization  problem,  derive  the  optimality  criteria,  and  verify  that  the 
numerical  results  match  the  derived  optimality  criteria.  The  problem  structure  is 
derived  in  Chapter  III,  followed  by  the  numerical  results  and  comparisons  in  Chapter 
IV. 
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III.  Problem  Definition  and  Assumptions 

This  chapter  documents  the  generic  dynamic  optimization  problem  as  well  as 
the  specific  constraints  used  throughout  this  research.  The  simpler  HCV  model  and 
mission  are  described  first.  This  is  used  to  build  up  to  the  more  complex  mission  using 
the  CAV  vehicle.  With  the  problem  and  assumptions  clearly  defined,  the  analysis  and 
results  follow  in  Chapter  IV. 

3.1  Generic  Problem  Statement 

The  overall  objective  of  this  research  is  to  create  an  optimal  trajectory,  for 
a  given  vehicle,  to  strike  a  final  target  in  minimum  time.  The  cost  or  objective 
function  defines  optimality,  which  is  mission  dependent  and  typically  defined  by  the 
user.  For  example,  the  user  may  want  to  arrive  in  minimum  time  or  maximum 
energy.  Additionally,  equality  and  inequality  constraints  are  required  to  define  the 
given  mission;  thus,  allowing  incorporation  of  waypoints  and  no-fly  zones.  Other 
path  limitations,  such  as  thermal  environment,  may  also  exist.  The  vehicle  itself 
is  represented  by  the  dynamics  or  equations  of  motion  and  any  control  limitations. 
These  generic  definitions  are  formulated  into  equations  and  used  to  solve  for  the 
optimal  control,  leading  to  the  optimal  trajectory.  Therefore,  the  first  step  is  to 
introduce  the  dynamic  optimization  nomenclature  used  throughout  this  research. 

3.1.1  Dynamic  Optimization.  The  continuous  time,  dynamic  optimization 
problem  with  open  final  time  is  the  particular  problem  of  interest.  The  open  final  time 
allows  time  to  be  in  the  cost  function,  which  is  the  intent  for  a  minimum  time  problem. 
The  advantage  to  the  continuous  time  formulation,  over  discrete  time  formulation, 
is  it  allows  the  direct  incorporation  of  the  differential  equations  of  motion.  The 
equations  of  motion  in  Equation  (1)  are  a  function  of  the  states,  x(f)  G  Mn,  the 
controls  u(t)  €  Mm,  and  the  independent  variable  f6  R;  where  n  is  the  number  of 
states  and  m  is  the  number  of  controls. 

*(t)  =  /(x(t),u(t),t)  (1) 


21 


In  Bolza  form  [12],  the  objective  or  cost  function,  J,  has  two  components.  The  first 
component  </>(x(fd), is  a  cost  defined  at  a  discrete  time  td,  or  the  sum  of  costs  at 
discrete  times,  0(x(trf1),  t^)  +  0(x(^2),  td2).  For  example,  for  final  time  tf  and  final 
state  x(ty)  =xj  the  first  term  is  0(xj,fy).  The  second  component  is  the  integrand 
cost  L ,  taken  from  the  initial  time  t0  to  final  time  tf.  The  combined  cost  function  is: 

J  =  0(x/,i/)  +  [  L(x(f),u(f),f)df  (2) 

Jto 

Additionally,  there  may  be  a  terminal  constraint  0  which  is  a  function  of  the  final 
state  and  the  final  time,  and  must  satisfy: 

^(x/,t/)=  0  (3) 

Certain  path  equality  constraints  C  may  also  exist  within  the  trajectory  [14:pg.99; 
44:pg.294],  Notice  here  that  the  equality  constraint  C  is  a  function  of  the  states  and 
the  controls: 

C(*(t),u(t),t)  =  0  (4) 

Inequality  constraints  are  also  important  to  this  research.  The  inequality  may  be  a 
function  of  either  the  states  or  the  controls  [14:pg.  108-1 18].  The  inequality  constraints 
only  appear  in  the  solution  during  the  times  they  are  active,  i.e.  equal  to  zero.  The 
state  inequality  constraint  is  S(x(£),£)  <  0.  When  the  control  constraint  is  active, 
the  control  inequality  constraint  in  Equation  (5)  is  a  special  case  of  Equation  (4), 
i.e.  C(u(t),t)  =  0.  Equation  (7)  illustrates  how  active  and  inactive  constraints  are 
handled.  The  technique  being  formulated  herein  requires  the  equality  constraint  to 
be  a  function  of  the  control;  therefore,  the  nomenclature  for  C  has  been  maintained 
whereas  the  distinction  is  made  for  S  since  it  is  not  a  function  of  the  control.  The 
separate  state  and  control  inequality  constraints  are: 

S'(x(t),t)<  0  C(u(f),f)<0  (5) 
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The  goal  of  optimization  is  to  determine  a  control  to  minimize  the  cost  while  satisfying 
the  dynamics  and  constraints.  In  order  to  ensure  the  constraints  are  satisfied  they 
are  adjoined  to  the  cost  function  using  Lagrange  multipliers;  u,  A,  and  fi.  Notice  that 
the  following  adjoined  cost  function  J  maintains  the  same  value  as  Equation  (2)  if  all 
of  the  constraints  are  satisfied.  To  simplify  notation,  the  dependence  of  each  variable 
is  omitted,  e.g.  0(x/,if)  is  simply  written  as  (j).  Thus  the  optimization  problem  is  to 
determine  the  control  u  that  minimizes  the  scalar  cost  J: 

min  J  —  (j)  +  uT^  +  f  [L  +  A T(/  —  x)  +  /iTC]  dt  (6) 

J  to 

To  activate  the  constraints  when  C  —  0,  the  applicable  components  of  the  multiplier 
/i  take  on  the  following  values: 


>0,  C  =  0, 

I-1-  < 

[  =  0,  C  <  0 

The  state  inequality  constraint  S  appears  to  be  missing;  however,  a  method  to  in¬ 
corporate  this  constraint  is  shown  later.  The  calculus  of  variation  dictates  that  the 
variations  of  the  adjoined  cost,  taken  with  respect  to  each  variable,  must  be  zero  in 
order  for  the  solution  to  possibly  be  a  minimum.  Thus,  the  variation  of  the  adjoined 
cost  function  must  be: 

5J  =  0  (8) 

Before  presenting  the  results,  some  intermediate  variables  are  defined.  The  Hamilto¬ 
nian  H  is  a  scalar  defined  as: 


H  =  L  +  \Tf  +  iiTC 


(9) 


and  $  is  the  scalar  defined  as: 


$  =  (f)  +  uTip 


(10) 
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The  time  derivative  of  <f>  is  defined  as: 


6  =  <1>X  x  +  4>t/  (11) 

The  following  summarizes  the  necessary  optimality  conditions  for  the  continuous  time, 
free  final  time  problem;  the  intermediate  steps  can  be  found  in  [12:pg.l59].  The 
dependence  notation  is  returned  here  to  clarify  the  conditions  that  are  functions  of 
time  versus  those  that  apply  only  at  the  final  time.  The  following  shortened  notation 
has  been  maintained;  x  is  a  simplification  for  x(t)  and  u  is  a  simplification  for  u(t). 
Lastly,  the  subscripts  represent  a  partial  derivative  with  respect  to  that  variable.  The 
costates,  A,  differential  equations  are: 

A T(t)  =  — iLx(x,  u,  t)  =  — Lx(x,  u,  t)  -  A T(f)/X(x,  u,  t)  -  imt (t)Cx(x.,  u,  t)  (12) 

The  final  condition  on  the  costates  are: 

=  'Mx/,*/)  =  0*(x/,i/)  +  vT^x{^f,tf)  (13) 

The  optimality  criterion: 

#u(x,  u,  t)  =  Lu(x,  u,  t)  +  A r/u(x,  u,  t)  +  LLT(t)Cu(x,  u,  t)  =  0  V  t  >  t0  (14) 
and  finally  the  transversality  condition: 


^(x/,^/)  +  L(X/,U/,^/)  =  0  (15) 

Notice  that  S  does  not  explicitly  appear  in  these  optimality  criteria;  however,  Sec¬ 
tion  3.1.3  covers  how  S  is  to  be  included. 

3.1.2  Discontinuous  Lagrange  Multipliers.  A  unique  situation  occurs  when 
there  are  interior-point  state  constraints,  i.e.  constraints  that  apply  at  a  single  un- 
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specified  time  within  the  trajectory.  This  brings  about  the  concept  of  discontinu¬ 
ous  Lagrange  multipliers.  These  discontinuities  will  also  be  called  “jumps”  in  the 
costates  [14:pg.  101] .  Assume  there  is  only  one  interior-point  constraint,  occurring  at 
the  unspecified  time  ti,  as  a  function  of  the  states  and  time: 

N(x(ti),ti)  =  Q  (16) 

The  time  ti  now  represents  the  final  time  of  the  previous  segment  and  the  initial 
time  of  the  next  segment;  thus,  let  fj"  represent  just  before  t\  and  tf  just  after  fy. 
Similar  to  the  final  condition  on  the  costates  in  Equation  (13),  there  is  now  a  jump 
to  accommodate  the  intermediate  condition  in  Equation  (16)  [14:pg.l03;  44:pg.278]: 

=  AT(iU  +  ;rJ^^  (17a) 

H(ti)  =  (17b) 

The  vector  of  multipliers  7 rn  must  be  solved  for  to  satisfy  Equations  (17a)  and  (17b). 
The  finite  number  of  interior-point  constraints  may  be  greater  than  1;  therefore, 
Equation  (17)  can  be  generalized  to  apply  to  any  one  of  the  jumps  simply  by  changing 
ti  to  the  applicable  time  of  the  interior-point  constraint,  e.g.  ti  or  tj. 

3.1.3  Path  Inequality  and  Equality  Constraints.  The  path  inequality  con¬ 
straint  is  introduced  in  Equation  (5)  as  a  function  of  the  state  and  time;  it  is  repeated 
here  for  convenience: 

S'(x(t), t)  <  0  (18) 

Since  the  control  does  not  explicitly  appear  in  Equation  (18),  the  time  derivative  of  S 
is  taken  until  the  control  does  appear.  This  is  demonstrated  in  the  Breakwell  problem 
in  Appendix  F.  If  it  takes  q  derivatives  for  the  control  to  appear,  this  time  derivative 
is  written  as  S™.  For  the  special  case  where  this  is  the  only  constraint,  then  C  =  S ^ 
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and  Equation  (9)  becomes: 


H  (x,  u,  t)  =  L(x,  u,  t)  +  A 7  /(x,  u,  t)  +  nT  u,  t)  (19) 

This  process  adds  interior-point  constraints  for  all  of  the  q—  1  derivatives  which  each 
must  equal  zero.  Thus,  the  interior-point  constraint  vector  in  Section  3.1.2  becomes: 

S(x(ii),ti) 

(x(f  i ) ,  1 1 ) 

This  interior-point  constraint  occurs  at  some  time  t\,  which  for  the  path  constraint, 
is  the  point  of  boundary  contact.  This  is  also  characterized  by  a  jump  in  the  costates 
occurring  when  the  constraint  becomes  active  at  time  t\ .  Since  is  a  function  of 
the  control,  it  may  be  possible  to  maintain  =  0;  this  represents  being  on  the 
constrained  arc.  Unlike  entering  the  constrained  arc,  there  is  no  jump  in  the  costates 
upon  exiting  the  constrained  arc.  Also,  since  is  a  function  of  the  states,  the 
partial  derivative  appears  in  the  propagation  of  the  costate  equations,  Equa¬ 
tion  (12).  The  optimal  control  continues  to  be  solved  from  Equation  (14).  If  the 
computed  optimal  control  would  lead  to  violating  the  constraint,  then  the  control 
to  remain  on  the  constraint  is  used  instead.  Therefore,  the  optimal  trajectory  may 
follow  the  constraint,  or  simply  touch  the  constraint,  based  on  whichever  minimizes 
the  Hamiltonian  without  violating  the  constraint. 

This  section  defined  the  continuous  time  dynamic  optimization  problem  with 
free  final  time.  Equality  and  inequality  constraints  were  addressed  as  a  function  of  the 
controls  or  the  states.  The  next  step  is  to  specify  the  particular  mission  requirements 
in  terms  of  constraints  to  be  incorporated  into  a  combined  dynamic  optimization 
problem  for  this  research.  The  following  sections  address  the  mission,  the  simpler 
HCV  model,  and  then  progress  to  the  more  complex  CAV  model. 


7V(x(ti),ti)  = 
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3.2  Mission  Assumptions 

The  overall  problem  has  mission  components  and  vehicle  components.  The  mis¬ 
sion  components  include  the  target,  waypoints,  and  no-fly  zones;  whereas  the  vehicle 
components  include  the  equations  of  motion,  control  limitations,  and  aerodynamic 
heating  limitations.  The  following  mission  assumptions  are  presented  first  to  begin 
to  scope  the  problem: 

1.  The  waypoints  are  specified  in  the  desired  sequence. 

2.  Waypoint  passage  must  be  directly  overhead. 

3.  Inner-loop  control  is  available.  Only  the  outer-loop  or  trajectory  generation  is 
addressed. 

4.  Waypoints  are  sufficiently  spaced  such  that  no  two  are  within  the  turn  radius 
of  the  vehicle. 

5.  Altitude  for  waypoint  passage  is  not  specified. 

6.  No-fly  zones  are  specified  as  circular  exclusion  zones  with  infinite  altitude. 

7.  No-fly  zones  must  not  be  violated. 

8.  Target  coordinates  and  final  altitude  are  specified. 

Item  1  ensures  the  user’s  mission  definition  is  not  altered.  Item  2  avoids  modifying 
the  solution  by  manipulating  the  penalty  of  a  near  miss.  Item  3  scopes  this  research 
effort  to  focus  on  the  primary  trajectory  optimization  objective.  Item  4  characterizes 
the  expected  performance  of  a  hypersonic  vehicle.  Item  5  attempts  to  avoid  over 
specifying  the  problem;  furthermore,  waypoint  altitude  does  not  enhance  the  mission 
objectives.  Item  6  simplifies  the  derivations;  however,  other  shapes  could  alteratively 
be  incorporated.  Item  7,  like  item  2,  avoids  modification  of  the  solution  by  manipu¬ 
lating  the  penalty  of  a  slight  penetration  of  a  no-fly  zone.  Item  8  ensures  the  solution 
does  not  place  the  vehicle  at  such  an  excessive  altitude  that  spike  heating  becomes 
a  factor  during  a  rapid  final  descent,  or  leave  the  vehicle  with  insufficient  altitude  to 
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perform  any  required  terminal  maneuvers.  The  following  sections  cover  the  vehicle 
specific  assumptions  and  individual  vehicle  descriptions. 

3.3  Hypersonic  Cruise  Vehicle  (HCV) 

3.3.1  Vehicle  Description.  As  mentioned  in  the  Introduction  the  HCV  is  a 
proposed  hypersonic  global  strike  vehicle.  For  this  study  its  thrust  provides  the  ability 
to  maintain  altitude  to  simplify  the  equations  of  motion  by  removing  the  variation 
in  altitude.  Thrust  is  also  used  to  maintain  constant  speed  or  maintain  constant 
deceleration.  For  this  research  thrust  is  not  a  control  variable,  it  is  merely  used 
conceptually  to  validate  the  constant  altitude,  speed,  and  deceleration  simplifications. 
The  maximum  available  thrust  will  dictate  the  achievable  flight  envelope,  i.e.  altitude 
and  velocity.  To  determine  the  validity  of  a  trajectory,  operational  limitations  are 
derived  from  a  similar  aircraft.  The  National  Aerospace  Plane  (NASP)  (UX-30)  [121] 
may  be  used  due  to  its  unclassified  data;  however,  the  SR-71  is  used  here  since  it  is 
a  proven  platform  and  the  data  is  readily  available  [101]. 

The  stated  objective  of  9,000  nautical  miles  (nmi)  in  2  hours  only  partially  de¬ 
fines  the  flight  envelope;  therefore,  some  addition  assumptions  are  made  and  Mach 
numbers  are  computed  for  comparison  to  other  aircraft.  A  representative  flight  alti¬ 
tude  of  100,000  ft  is  chosen  since  it  was  the  approximate  test  altitude  for  the  unmanned 
hypersonic  Hyper-X  (X-43A)  test  vehicle  [98].  The  range  and  time  objectives  equate 
to  2.315  km/s,  which  at  flight  level  (FL1000)  is  Mach  7.66.  In  the  Operating  Lim¬ 
itations  section  of  the  SR-71  Flight  Manual  [101],  equivalent  airspeed,  expressed  in 
knots  equivalent  airspeed  (KEAS),  is  used  to  specify  a  minimum  acceptable  airspeed 
over  a  range  of  altitudes  and  Mach  numbers.  This  minimum  is  used  to  verify  suffi¬ 
cient  airspeed  to  sustain  level  flight  or  engine  performance.  This  minimum  acceptable 
airspeed  is  Vrnm  =  310  KEAS,  which  equates  to  Mach  4.47  at  FL1000.  If  a  trajec¬ 
tory  reaches  the  target  at  an  airspeed  less  than  Vmin  it  is  considered  invalid  since  the 
constant  altitude  assumption  has  been  violated.  For  example,  a  flight  starting  at  an 
initial  airspeed  of  2.315  km/s  at  FL1000,  as  described  above,  with  a  deceleration  of  0.1 


m/s2  over  2  hours  has  a  final  airspeed  of  366  KEAS  or  Mach  5.28,  thus  the  level  flight 
assumption  remains  valid.  A  deceleration  profile  is  chosen  to  demonstrate  the  effects 
of  non-constant  velocity  which  leads  to  non-constant  turn  radii.  Also,  deceleration 
will  be  a  component  of  the  more  complex  CAV  model.  The  bank  angle  a  is  limited 
to  ±20°  due  to  wing  loading,  stability,  and/or  inlet  airflow  tolerances.  A  constant 
altitude  is  maintained  by  ensuring  the  vertical  component  of  lift  equals  the  vehicle 
weight  i.e.  Lift  =  mg/  cos  a.  The  thrust  used  to  counteract  the  drag  associated  with 
the  additional  lift  due  to  a  non-zero  bank  angle  is  Thrust  =  m  (a  +  g/(cos  ctCl/Cd)). 
These  derived  parameters  are  used  as  a  representative  HCV  for  demonstrating  results. 

3.3.2  Mission  Profile.  For  the  analysis,  a  particular  mission  is  specified 
in  Table  1  by  its  initial  location,  intermediate  waypoints,  no-fly  zone(s),  and  final 
destination  or  target.  The  waypoints  represent  specified  mission  critical  locations 
that  must  be  flown  over  precisely,  either  for  reconnaissance  or  for  multiple  payload 
deliveries.  It’s  arguable  that  a  vehicle  may  be  limited  to  level  flight  for  imaging  or 
payload  deployment;  however,  due  to  the  long  durations  of  the  turns,  a  short  level 
segment  within  the  turn  is  approximated  here  as  a  continuous  turn.  Therefore,  in 
addition  to  no  heading  angle  constraint,  there  is  also  no  bank  angle  constraint  at  the 
waypoints.  The  no-fly  zone  coordinates  are  the  center  of  the  keep  out  circle  with 
radius  specified. 


Table  1:  HCV  Mission  Description 


Descriptor  Latitude  Longitude  Radius 


Initial 

N 

28° 

34' 

W 

80° 

38' 

Waypoint  1 

N 

52° 

11' 

W 

43° 

46' 

Waypoint  2 

N 

17° 

01' 

W 

20° 

33' 

No-Fly  Zone 

N 

22° 

39' 

E 

11° 

03'  840  nmi 

Target 

N 

33° 

46' 

E 

29° 

34' 

fio=100,000  ft=30.48  km,  V0=2.315  km/s,  0o=O° 


29 


The  velocity  is  V ;  and  the  heading  angle  is  9,  measured  positive  going  counter¬ 
clockwise  from  the  ;r-axis. 

3.3.3  HCV  Assumptions.  The  assumptions  made  for  the  HCV  simplify  the 
analytic  derivation  of  the  optimality  conditions.  This  simplification  makes  it  easier  to 
follow  the  derivation,  and  helps  add  intuition  into  the  final  results.  The  assumptions 
are: 

1.  Altitude  is  maintained  regardless  of  velocity. 

2.  Deceleration  is  either  zero  or  constant. 

The  first  assumption  enforces  constant  altitude;  however,  post  run  analysis  assesses 
the  validity  of  the  flight  profile  based  on  an  acceptable  minimum  velocity.  The  second 
assumption  decouples  the  bank  angle  and  changes  in  velocity.  Realistically,  an  increase 
in  bank  angle  would  decrease  the  vertical  component  of  lift  and  the  angle-of-attack 
would  have  to  increase  to  provide  enough  lift  to  maintain  altitude.  The  increase  in 
angle-of-attack  would  then  result  in  an  increase  in  drag  resulting  in  a  decrease  in 
velocity.  Again,  the  intent  of  this  simpler  model  is  to  gain  insight  into  the  dynamic 
optimization  problem,  not  to  produce  a  high  fidelity  ffCV  model. 

3.3.4  Cost  Function.  The  objective  is  to  minimize  flight  time.  The  cost  can 
therefore  be  expressed  as  tj  or  as  f*f  ldi.  The  former  seems  to  produce  slightly  faster 
results;  possible  because  the  integral  of  1  internally  requires  an  additional  state.  The 
results  should  obviously  be  the  same  with  either  expression,  so  for  this  research  the 
cost  function  is: 

J  =  tf  (21) 

3.3.5  2-D  Equations  of  Motion.  The  2-D  equations  of  motion  confine  the 
flight  profile  to  a  constant  altitude,  horizontal  plane.  Also,  the  control  u  is  a  function 
of  the  bank  angle  normalized  by  the  tangent  of  the  maximum  bank  angle  (crmax),  u  = 
tan  a  j  tan  amax.  The  purpose  for  the  normalization  is  to  easily  identify  the  minimum 
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and  maximum,  i.e.  =|=1.  The  units  for  the  states  and  time  are  nondimensional  in  order 
to  facilitate  numerical  stability.  The  conversion  from  dimensional  to  nondimensional 
is  in  Appendix  B  as  Equation  (B.17).  The  2-D  HCV  equations  of  motion  are: 


x 

y 

6 

V 

Where  a  is  a  constant.  Now  dehne 

x  =  x  y 


V  cos  9 

(22a) 

V  sin  9 

(22b) 

tan  (7 max 

v  u 

(22c) 

a 

(22d) 

1  T 

V  =/(x,u,t) 

(23) 

3.3.6  Control  Constraints.  For  this  model  of  the  HCV  the  only  control  is 
the  normalized  bank  angle  u.  As  described  in  Section  3.3.1,  the  maximum  bank  angle 
is  ±20°.  In  Section  3.3.5  the  control  is  normalized  to  produce  an  acceptable  range  for 
u,  i.e.  —1  <  u  <  1.  To  dehne  an  equality  constraint  as  a  function  of  the  control,  as 
described  in  Section  3.1.1,  the  control  constraint  is: 


u  —  1 
— u  —  1 


<  0 


(24) 


3.3.7  Terminal  Constraints.  The  terminal  constraint  for  the  HCV  are 
merely  the  specified  final  coordinate  [xf,  i/f],  i.e.  the  target.  Altitude  does  not  change 
so  it  does  not  need  to  be  included: 


V>(x(t/),t/) 


x(tf)  -  Xf 

y(tf )  -  yf 


(25) 
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3.3.8  Path  Inequality  Constraints.  The  no-fly  zone  assumption  in  Section  3.2 
states  that  the  no-fly  zones  are  circular  exclusion  zones.  Since  there  can  be  multiple 
110-fly  zones  each  one  is  indexed  j.  Each  no-fly  zone  is  characterized  by  the  coordinate 
of  its  center,  [xCj,  yCj ],  and  its  radius,  Rj.  Thus,  to  remain  outside  the  no-fly  zone 
the  current  distance  from  the  center  of  the  no-fly  zone  must  be  equal  to  or  greater 
than  the  110-fly  zone  radius.  Let  the  displacement  in  the  ^-direction  and  ^-direction 
be,  A Xj  =  x(t)  —  xCj  and  A yj  =  y(t)  —  yCj,  respectively.  In  equation  form,  consistent 
with  Section  3.1.3,  each  no-fly  zone  inequality  constraint  is: 

S(xW.f)  =  1  {R]  _  ^  <  o  (26) 

J  —  ^  •••)  Jend  ^ 

The  value  jend  is  the  user  specified  finite  number  of  no-fly  zones.  This  constraint  is 
in  effect  for  all  time  t:  however,  the  next  section  shows  there  is  also  an  unspecified 
discrete  time  tj  associated  with  each  no-fly  zone. 


3.3.9  Interior- Point  Constraints.  The  waypoints  and  no-fly  zones  both 
create  interior-point^  constraints.  The  waypoints  are  already  the  definition  of  an 
interior-point  constraint;  that  is  a  specified  function  of  the  states  at  a  discrete  un¬ 
specified  time  that  must  equal  zero.  Since  there  may  be  multiple  waypoints,  each  one 
is  indexed  i,  thus  the  individual  waypoint  passage  occurs  at  time  p.  The  constraint 
to  fly  over  the  waypoint  at  time  ti  implies  the  vehicle  coordinates,  [x(ti),  y(ti)],  are 
equal  to  the  specified  waypoint  coordinates,  [xi,  yi\.  The  interior-point  constraints, 
for  the  finite  number  of  waypoints  i  to  iend,  are  thus: 


N(x.(ti),ti) 


i  1,2,...,  iend 


X(ti)  ~  Xi 

y(U )  -  Vi 


(27) 


'  An  interior-point  within  this  research  corresponds  to  a  time  along  the  trajectory  between  the 
initial  and  final  time,  thus  creating  a  multi-point  boundary  value  problem.  Hence,  as  referenced 
herein,  an  interior-point  constraint  is  not  associated  with  the  interior  point  method  used  within  the 
field  of  linear  programming. 
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Section  3.1.3  demonstrates  how  the  110-fly  zone  inequality  constraint  can  produce 
interior-point  equality  constraints.  Specifically,  the  time  derivative  of  the  inequality 
constraint,  Equation  (26),  must  be  taken  until  a  control  in  u  appears.  For  the  HCV, 
u  appears  in  the  second  derivative;  therefore,  the  interior-point  constraint  associated 
with  each  no-fly  zone  consists  of  the  zero  derivative  S  and  the  first  time  derivative 
.  These  are  tangency  requirements,  meaning  the  boundary  must  be  approached 
at  a  tangent  or  a  violation  will  occur.  This  takes  the  same  form  as  the  waypoint 
constraints  N  but  M  will  be  used  for  distinction: 

3  1)2,...,  Jend 

Also,  the  multipliers  in  Equation  (17)  are  renamed  7rm: 


AT(t-) 

.T,  T  DM 

=  *  (t+)+^9x(t.) 

(29a) 

m) 

TT,  Td M 

=  H{f>)  ^gt. 

(29b) 

This  equality  constraint  S  differs  from  the  inequality  constraints  in  Equation  (26) 
since  this  S  only  applies  at  the  time  of  no-fly  zone  boundary  contact  tr  The  next 
section  handles  the  second  derivative  of  the  no-fly  zone  constraint  which  must  be  zero 
to  remain  on  the  no-fly  zone. 

3.3.10  Path  Equality  Constraints.  As  derived  in  Section  3.1.1,  a  state  equal¬ 
ity  constraint  is  adjoined  to  the  Hamiltonian,  Equation  (19).  For  this  problem  a  state 
equality  constraint,  generically  expressed  as  C(x,u,t)  in  Section  3.1.1,  is  required  to 
remain  on  the  no-fly  zone  constraint  arc  until  returning  to  the  optimal  bang-level- 
bang  control  law;  specifically,  the  second  derivative  of  the  no-fly  zone  must  remain 


S(x(tj),tj) 

S(1)(x(tj),tj) 


I  (R*  -  Ax*  -  Ay*) 

-A XjV  cos  6  —  Ay3 V  sin  0 


=  0  (28) 
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zero  until  departing  the  no-fly  zone  boundary: 


(7(x,  u,  t)  =  (x,  u,  t) 

=  —V2  +  u(Axj  sin  9  —  A yj  cos  9)  tan  omax  +  (—A Xj  cos  9  —  A yj  sin  6) a 
=  0  (30) 

The  partial  derivative  is  used  to  propagate  the  costates. 

3.3.11  HCV  Summary.  The  cost,  dynamics,  control  limitations,  waypoints, 
and  no-fly  zones  have  all  been  expressed  in  terms  of  a  dynamic  optimization  problem. 
Three  types  of  state  equality  constraints  have  been  shown;  interior-point  constraints, 
path  equality  constraints  or  tangency  requirements,  and  maintaining  boundary  con¬ 
tact.  The  principal  difference  is  interior-point  constraints  and  tangency  requirements 
occur  at  a  single  time,  whereas  the  adjoined  q  time  derivative  path  constraints  apply 
the  entire  time  the  particular  constraint  is  active,  i.e.  equal  to  zero.  This  distinction 
dictates  the  form  of  the  solution  as  either  a  jump  in  the  costates,  as  in  Equation  (17), 
or  a  propagation  of  the  costates,  as  in  Equation  (12).  The  analytical  solution  to  this 
simpler  HCV  problem  is  found  in  Chapter  IV,  and  the  more  complex  CAV  problem 
setup  is  presented  in  the  next  section. 

3.4  Common  Aero  Vehicle  (CAV) 

The  CAV  is  an  operationally  representative  vehicle  for  this  research.  It  rep¬ 
resents  a  hypersonic  reentry  vehicle  with  crossrange  capability  necessary  to  acquire 
waypoints  and  avoid  no-fly  zones.  It  has  no  thrust;  therefore,  the  reentry  profile 
must  avoid  exceeding  a  specified  limit  of  the  heating  rate  at  the  stagnation  point. 
The  following  describes  the  vehicle  model  and  expresses  the  CAV  mission  in  terms  of 
dynamic  optimization  constraints. 

3-4-1  Vehicle  Description.  Reference  [69]  describes  a  low-lift  CAV  and  a 
high-lift  CAV.  The  high-lift  CAV,  CAV-H,  is  modeled  here  to  extend  the  crossrange 
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capability.  The  full  aerodynamic  database  and  characteristic  parameters  are  presented 
in  Appendix  A.  The  control  and  heating  limitations  are  realistic,  but  are  chosen  to 
force  boundary  contact  in  order  to  verify  different  forms  of  the  solution. 

3-4-2  Mission  Profile.  A  target  nearly  halfway  around  the  globe  is  used  to 
represent  a  global  strike  mission.  A  waypoint  in  the  middle  of  the  Atlantic  represents 
a  telemetry  station  that  can  validate  the  vehicle  navigation  and  control  to  terminate  a 
rogue  vehicle  if  necessary.  The  radius  of  the  first  no-fly  zone  is  intentionally  chosen  to 
be  much  smaller  than  the  turn  capability  of  the  CAV,  thus  forcing  the  solution  to  only 
be  able  to  contact  the  boundary  at  a  single  point.  The  next  waypoint  represents  a 
secondary  mission  target  for  reconnaissance  or  payload  delivery.  The  last  no-fly  zone 
has  a  large  enough  radius  to  allow  the  boundary  to  be  followed  if  optimal,  and  is  large 
enough  to  force  full  control  authority  in  order  to  clear  the  no-fly  zone  and  still  make 
it  to  the  target.  The  heating  constraint  limitation  is  set  low  enough  such  that  the 
unconstrained  heating  problem  exceeds  the  value,  thus  the  constrained  solution  must 
incorporate  the  limitation.  The  initial  conditions  are  approximately  those  presented 
in  [69] .  The  specified  final  altitude  is  the  approximate  altitude  that  the  Space  Shuttle 
uses  to  complete  its  entry  guidance  phase  [39].  Table  2  provides  the  mission  objectives: 


Table  2:  CAV  Mission  Description 


Descriptor 

Latitude 

Longitude 

Radius 

Initial 

N  28°  35.286' 

W  80°  40.194' 

Waypoint  1 

N  34°  2.810' 

W  27°  18.430' 

No-Fly  Zone  1 

N  20°  15.513' 

W  3°  27.588' 

960  nmi 

Waypoint  2 

N  33°  13.298' 

E  41°  41.266' 

No-Fly  Zone  2 

N  55°  43.849' 

E  58°  33.688' 

1500  nmi 

Target 

N  31°  36.653' 

E  65°  42.016' 

h0= 122  km,  V0- 

=24,000  ft/s=7.3152  krn/s,  7o=- 

■1.5°  6*0=4° 
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Altitude  is  h.  The  flight  path  angle  7,  as  depicted  in  Figure  12  [41],  is  the  angle 
with  respect  to  the  local  horizontal,  measured  positive  away  from  the  Earth. 

1/ 


Figure  12:  Flight  Path  Angle  Defined  [41] 

3-4-3  CAV  Assumptions.  The  following  assumptions  simplify  the  3-D  equa¬ 
tions  of  motion  and  maintain  the  x,  y,  and  6  states  presented  for  the  2-D  case.  These 
simplifications  do  not  indicate  limitations  to  this  solution  process,  they  are  intended 
to  make  the  derived  analysis  more  intuitive  and  easier  to  follow.  The  CAV  3-D  model 
assumptions  are: 

1.  Flat,  non-rotating  Earth 

2.  Gravity  is  constant 

3.  Flight  path  angle  is  small 

4.  Drag  is  the  dominate  deceleration  term 

5.  Coefficient  of  Lift  (Cl)  and  Coefficient  of  Drag  ( Cd )  are  only  a  function  of 
angle-of-attack 

6.  Control  is  bank  angle  and  angle-of-attack,  both  limited 

7.  Atmospheric  density  is  modeled  as  a  simple  exponential 

8.  bleat  flux  at  the  stagnation  point  is  limited 

9.  G-loading  and  total  heat  are  not  addressed 
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Ignoring  the  rotation  of  the  Earth  is  a  common  assumption  seen  in  [41;  97;  109;  110]. 
The  flat  Earth  approximation  is  most  applicable  to  lifting  bodies  with  a  low  flight  path 
angle  as  seen  in  [114:pg.248;  87;  73:pg.204] .  The  flat  Earth  model  is  used  as  a  proof 
of  concept  and  to  provide  commonality  with  the  2-D  derivation;  however,  a  spherical 
Earth  model  would  be  more  appropriate  for  computing  operationally  representative 
results.  The  small  flight  path  angle  leads  to  the  small  angle  approximations.  For 
hypersonic  vehicles,  Cl  and  Cn  are  often  assumed  independent  of  Mach  Number,  i.e. 
simply  a  function  of  angle-of-attack  [111].  The  estimated  fit  to  the  CAV  aerodynamic 
data  is  in  Appendix  A. 

3-4-4  Cost  Function.  The  minimum  time  objective  function  from  the  HCV 
in  Section  3.3.4  remains  the  cost  function  for  the  CAV: 

J  =  tf  (31) 

3-4-5  3-D  Equations  of  Motion.  The  equations  of  motion  represent  a  flat 
Earth  model  with  zero  Earth  rotation.  A  small  angle  approximation  for  the  flight 
path  angle  implies  cosy  ~  1  and  sin 7  rs  7.  It  is  also  assumed  that  the  drag  term  is 
the  dominate  term  in  the  V  equation,  i.e.  the  component  of  gravity  in  the  velocity 
direction  is  negligible  compared  to  the  aerodynamic  drag.  This  last  assumption  is 
consistent  with  the  small  angle  approximation  for  7.  The  nondimensional  equations 
of  motion  originating  from  [114:pg.249;  87]  and  derived  in  Appendix  B  are: 


X  = 

V  cos  6 

(32a) 

y  = 

V  sin  6 

(32b) 

h  = 

V7 

(32c) 

v  = 

By2e-proh  +  c2) 

(32d) 

2  E* 

7  = 

BVe-prohctcosa  -  ^  +  V 

(32e) 

6  = 

BVe-^ce  sin  a 

(32f) 

37 


The  constant  B  is  a  function  of  the  vehicle  parameters,  atmospheric  constant,  and 
initial  conditions.  The  constant  E*  is  the  maximum  lift-to-drag  ratio.  The  controls 
are  bank  angle,  cr,  and  the  fraction  of  C*L)  q.  The  coefficient  C *L  is  the  coefficient 
of  lift  that  produces  the  maximum  lift-to-drag  ratio.  More  details  are  contained  in 
Appendix  B. 


3-4-6  Control  Constraints.  The  controls  are  bank  angle,  a,  and  normalized 
coefficient  of  lift,  q.  The  bank  angle  is  limited  to  ±60°  to  represent  stability  lim¬ 
itations,  and  to  force  the  control  to  reach  the  imposed  limit  for  this  research.  The 
maximum  limitation  of  q  =  2  is  derived  from  the  aerodynamic  data  in  Appendix  A. 
The  minimum  q  =  0  maintains  proper  orientation  for  the  thermal  protection  system. 
Thus  the  control  limitations  are: 


C{u,t) 


<J  —  7r/3 
—a  —  7r/3 
Cl-  2 


<  0 


-ci 


(33) 


3-4-7  Terminal  Constraints.  The  terminal  constraints  for  the  CAV  are  the 
target  coordinates  and  the  specified  final  altitude  hf. 


x(tf)  -  xf 

y(tf)  -  vj 

h(tf )  -  hf 


(34) 


3-4-8  Path  Inequality  Constraints.  The  path  inequality  constraints  for  the 
CAV  are  the  no-fly  zones  and  the  maximum  stagnation  point  heat  rate.  The  no-fly 
zone  discussed  in  Section  3.3.8  is  repeated  here  for  convenience: 


j  1)2,...,  Jend 


-  [R]  -  Axj  -  Ay])  <  0 


(35) 
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The  function  for  the  heating  constraint  is  found  in  [41;  109],  which  has  its  origin 
from  [15],  and  is  also  discussed  in  [34],  The  dimensional  form  of  the  equation  for 
Qsdim  IS. 


Qsdim 


k 


dim 


p 

Psl  , 


y/2  /  Yam  V 
/  VVS»ro/ 


(36) 


The  dimensions  are  based  on  the  converted  units  of  the  constant  kdim,  which  is  17,000 
Btu  ft-3/2s-1  [15],  and  vehicle  dependent  radius  of  the  nose,  rnose.  The  dimensional 
maximum  heat  flux,  qsmax,  is  used  to  nondimensionalize  Equation  (36).  The  nondi- 
mensional  altitude  h  and  velocity  V  are  defined  in  Appendix  B.  The  nondimensional 
heating  rate  at  stagnation  point,  qs,  in  terms  of  the  nondimensional  states  h  and  V 
is: 


qs  =  Ke-Pr0h/2V3  (37) 

where  K  =  kdim^~^ro~Rs>^2 / (\Anose  Qs  max )  making  the  maximum  heat  flux  qs  =  1. 
The  heating  path  inequality  constraint  will  have  the  same  form  as  S  in  Equation  (35); 
however,  Q  is  used  for  distinction.  There  is  no  minimum  heating  limitation;  however, 
from  Equation  (37)  qs  is  always  >  0.  Thus  the  heating  inequality  constraint  is: 


Q(x(t),  t)  =  [qs  -  1]  =  [Ke~/3r°hl2V3  -  l]  <  0  (38) 


This  Q  is  not  related  to  total  heating  which  is  often  represented  by  this  symbol.  For 
this  research  the  total  heating  is  considered  a  vehicle  design  parameter  and  not  a 
trajectory  optimization  parameter. 


3-4-9  Interior-Point  Constraints.  Since  the  flat  Earth  model  maintains  the 
states  x  and  y,  the  waypoint  interior-point  constraint  from  Equation  (27)  still  applies: 


N(xL(ti),ti) 
i  1,2,...,  lend 


x{ti)  -  Xi 

y(U )  ~  Vi 


(39) 
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Similarly,  the  derived  no-fly  zone  interior-point  constraint,  also  from  Section  27,  is 

still: 


j  =  1)2,...,  Jend 

It  can  be  seen  that  the  heating  constraint,  Equation  (38),  is  not  a  function  of  either 
control.  Therefore,  as  described  in  Section  3.1.3,  time  derivatives  must  be  taken 
until  a  control  explicitly  appears.  This  only  requires  the  first  time  derivative,  thus 
the  interior-point  constraint  derived  from  the  heating  constraint  Q  up  to  the  q  —  1 
derivative  is  simply: 


1  (R2  -  Ax2  -  Ay |) 
-AxjV  cos  9  —  A yjV  sin  9 


=  0  (40) 


k 


Q(x(4),4) 
1,2,...,  ke nd 


[qs(h)  -  1]  =  [Ke~^h^/2V{tkf  -  1]  =  0 


(41) 


This  equation  only  requires  satisfaction  upon  contact  with  the  heating  constraint  at 
time  tk,  where  k  represents  the  number  of  contacts  with  the  heating  constraint. 


3-4- 1 0  Path  Equality  Constraints.  The  path  equality  constraints  apply  while 
a  constraint  is  active,  i.e.  while  the  inequality  constraint  equals  zero.  The  control  to 
remain  on  the  constraint  arc  is  used  if  the  optimal  control  cannot  be  maintained  due 
to  constraint  violation.  Therefore,  if  on  the  no-fly  zone  constraint,  the  second  time 
derivative  of  the  inequality  constraint  must  equal  zero: 


(7(x,  u,  t)  =  (x,  u,  t) 

=  —V2  +  (A Xj  sin  9  —  Ayj  cos  9)  6  +  (—A Xj  cos  9  —  Ayj  sin  9)  V 
=  0  (42) 


The  form  of  this  constraint  differs  from  the  HCV  problem  since  the  CAV  equations 
of  motion  are  different.  The  term  in  front  of  V  is  a  multiple  of  thus  will  be 
zero;  however,  it  is  retained  in  order  to  take  the  partial  derivative  with  respect  to 
x  to  propagate  the  costates.  Next,  while  on  the  heating  constraint,  the  first  time 
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derivative  of  the  heating  inequality  constraint  must  equal  zero: 


Q(1)(x,t)  =  +  —  e~^\l  +  cj)  =  0  (43) 

Once  again,  contact  with  a  constraint  does  not  dictate  remaining  on  the  constraint. 
A  boundary  arc  is  only  traversed  if  the  optimal  control  would  result  in  violation  of 
the  constraint.  This  will  be  analyzed  in  the  following  Analysis  chapter. 

3-4- 11  CAV  Summary.  The  CAV  model  increases  the  model  complexity 
with  additional  degrees-of-freedom;  namely,  altitude  and  flight  path  angle.  The  vari¬ 
able  altitude  leads  to  a  stagnation  heating  constraint.  While  the  waypoint  and  no-fly 
zone  constraints  have  the  same  form,  their  analytical  solutions  are  different  than  the 
HCV.  The  differences  arise  from  the  derived  constraints  being  functions  of  the  state 
derivatives;  therefore,  since  the  equations  of  motion  have  changed,  the  derived  con¬ 
straint  equations  change.  This  will  be  covered  further  in  the  Analysis  chapter  to 
follow. 

3.5  Numerical  Methods  Implementation  Techniques 

The  next  chapter  will  show  that  the  simpler  HCV  problem  has  an  analyti¬ 
cal  solution;  contrarily,  the  complexity  of  the  CAV  requires  numerical  techniques  to 
efficiently  compute  the  solution.  A  numerical  technique  typically  represent  the  con¬ 
tinuous  problem  via  some  form  of  discretization.  The  discretization  method,  i.e.  the 
spacing  of  the  discretized  points,  is  going  to  be  based  on  the  dynamics.  Thus,  the 
discretization  may  need  to  be  different  for  different  portions  of  the  modeled  dynam¬ 
ics.  For  example,  Runge-Kutta  integration  is  often  run  with  variable  step  size  to 
allow  for  an  increase  in  the  number  of  time  steps  during  rapid  state  changes.  This 
methodology  remains  true  when  attempting  to  model  the  original  states  plus  the 
added  dimensionality  of  the  costates;  therefore,  an  increased  number  of  time  steps 
is  required  for  rapid  changes  in  the  states  as  well  as  rapid  changes  in  the  costates. 
Section  3.1.2  discusses  discontinuous  Lagrange  multipliers,  i.e.  jumps  in  the  costates. 
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Thus,  in  order  to  properly  model  the  system,  any  numerical  technique  must  be  able 
to  model  these  discontinuities  accurately.  This  modeling  can  be  achieved  by  creating 
a  dense  set  of  time  steps,  or  by  actually  allowing  for  a  discontinuity.  The  dense  set 
of  time  steps  is  used  when  the  numerical  method  forces  continuity,  thus  the  dense  set 
of  continuous  points  replicates  a  discontinuity.  Alternatively,  a  model  can  allow  for 
discontinuities  for  certain  parameters  while  maintaining  continuity  for  others.  For  ex¬ 
ample,  the  controls  or  costates  may  have  permissible  discontinuities  whereas  the  states 
may  be  forced  to  be  continuous.  These  difference  in  discretization  are  characteristic 
of  the  software  packages  DIDO  and  GPOCS.  DIDO  increases  the  density  of  the  time 
steps,  called  nodes,  by  specifying  intermediate  “knots”,  while  GPOCS  allows  the  user 
to  break  the  problem  up  into  phases.  The  phased  approach  of  GPOCS  allows  for 
a  multi-point  boundary  value  problem  (MPBVP)  by  specifying  terminal  constraints 
at  the  end  of  each  phase,  i.e.  interior-point  constraints.  Additionally,  each  phase  is 
discretized  individually;  therefore,  the  Chebyshev  spacing  described  in  Appendix  E 
creates  the  desired  situation  of  higher  node  density  at  the  beginning  and  end  of  each 
phase.  The  choice  of  phase  breakpoints  dictates  where  desired  higher  accuracy  oc¬ 
curs  or  where  discontinues  are  allowed.  This  Chebyshev  spacing  can  have  the  adverse 
affect  of  leaving  sparse  spacing  in  the  middle  of  the  phase,  even  with  an  increase  in 
the  specified  number  of  nodes.  If  this  middle  region  requires  higher  accuracy,  it  may 
be  necessary  to  force  another  phase  change  in  this  region.  In  summary,  the  CAV 
numerical  trajectory  solution  is  broken  up  into  phases  for  three  purposes;  to  specify 
interior-point  constraints,  to  allow  for  discontinuities,  and/or  to  increase  the  accuracy 
at  a  region  of  high  state  or  costate  change. 

One  of  the  challenges  of  a  phased  trajectory  solution  or  approach  is  determining 
the  required  number  of  phases.  To  maximize  accuracy,  every  jump  in  the  costates 
should  occur  at  a  phase  breakpoint  in  order  to  allow  for  such  discontinuities.  The 
number  of  jumps  is  dictated  by  the  number  of  times  the  path  constraints  are  contacted; 
unfortunately,  this  number  is  typically  part  of  the  solution  and  is  not  known  a  priori. 
It  is  possible  to  run  a  lower  number  of  phases,  then  add  phase  breaks  at  each  path 
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constraint  contact;  however,  it  may  also  be  possible  to  allow  the  numerical  solution  to 
model  these  discontinuities  with  continuous  costates,  i.e.  without  a  phase  break.  For 
simplified  end  user  implementation,  it  is  therefore  desirable  to  minimize  the  number  of 
phases,  i.e.  eliminate  the  need  to  add  more  breakpoints.  The  ability  to  still  converge 
on  the  optimal  solution  even  with  a  minimal  number  of  phases  is  therefore  one  of 
the  topics  of  this  research.  Section  4.10  provides  a  comparison  between  a  trajectory 
broken  up  into  a  higher  number  of  phases,  i.e.  increased  accuracy,  and  a  trajectory 
broken  up  into  a  minimal  number  of  phases,  i.e.  more  desirable  user  implementation. 

3. 6  Solution  Comments 

This  section  provided  vehicle  and  mission  descriptions  for  both  the  ffCV  and 
the  CAV.  It  also  outlined  some  of  the  required  analytical  and  numerical  techniques 
to  finally  compute  a  solution  for  the  optimal  reentry  trajectory.  The  following  chap¬ 
ter  provides  the  detailed  analysis  and  numerical  results  for  both  vehicles  and  their 
respective  missions.  The  results  are  not  merely  presented,  but  are  broken  down  to 
ensure  they  match  the  expected  analytical  derivations.  User  strategies  for  repeatable 
solution  convergence  are  also  described.  The  outcome  is  a  verified  optimal  solution 
with  detailed  user  implementation  techniques. 
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IV.  Analysis  and  Results 


The  culmination  of  this  research  effort  lies  in  the  analysis  and  results  presented 
in  this  chapter.  Chapter  III  developed  the  fundamental  theory;  however,  the  problem 
specific  solutions  remain  to  be  solved.  For  each  vehicle  and  mission,  the  goal  is  to 
analytically  compute  the  optimal  solution  to  the  maximum  extent  possible.  For  the 
two-dimensional  constant  speed  HCV,  a  complete  geometric  solution  is  found  requir¬ 
ing  no  integration.  When  deceleration  is  added  to  the  HCV  problem,  a  solution  is 
determined  using  an  indirect  method  by  satisfying  an  analytically  derived  set  of  op¬ 
timality  criteria  [48].  Lastly,  a  proposed  direct  numerical  solution  technique  is  used 
to  solve  the  optimal  three-dimensional  CAV  reentry  trajectory  problem.  A  two-step 
process  is  implemented  to  validate  these  numerical  results.  First,  the  solution  to  the 
HCV  problem  is  recomputed  using  this  new  numerical  solution  technique  and  the 
results  are  compared  to  those  previously  computed  via  the  geometric  and  analytical 
approach.  Second,  the  candidate  numerical  result  is  tested  to  verify  it  matches  the 
derived  analytical  form  of  the  solution.  Once  this  exhaustive  process  of  verifying  the 
optimality  of  the  results  is  complete,  a  minimal  set  of  user  representative  steps  is  out¬ 
lined  for  future  implementation.  This  will  fulfil  the  research  objective  of  determining 
a  solution  process  capable  of  rapidly  converging  on  the  optimal  reentry  trajectory 
solution  that  satisfies  the  vehicle  dynamic,  heating  limitation,  and  the  added  mili¬ 
tary  constraints  of  waypoints  and  no-fly  zones;  all  within  the  confines  of  the  available 
control  authority. 

4- 1  2-D  Baseline 

In  order  to  access  the  improvement  achievable  by  optimizing  a  trajectory,  a 
lion-optimal  baseline  trajectory  is  computed  for  comparison.  This  baseline  trajec¬ 
tory  requires  very  little  computation  since  it’s  a  simple  steer  and  point  approach  to 
navigation.  The  baseline  controller  is  the  optimal  bang-level-bang  controller  for  each 
segment  separately  [23;  26;  88];  however,  this  segment  by  segment  technique  does 
not  ensure  optimality  over  the  entire  trajectory.  Specifically,  due  to  the  lion-optimal 
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incoming  heading  to  each  waypoint,  it  is  likely  to  overshoot  thus  extending  the  trajec¬ 
tory.  The  amount  of  time  to  be  gained  by  optimizing  over  this  baseline  is  proportional 
to  the  amount  of  time  the  vehicle  is  in  a  turn  compared  to  the  total  range.  For  a 
subsonic  aircraft,  with  relatively  small  turn  radii,  the  ratio  of  time  within  a  turn  is 
low  and  this  optimization  may  produce  a  negligible  time  savings.  However,  the  HCV 
has  huge  turn  radii  in  comparison  to  its  designed  range  and  thus  the  time  savings  can 
be  significant.  This  minimum  time  solution  will  likely  also  benefit  a  fighter  aircraft 
in  confined  terrain  or  a  reconnaissance  UAV  with  a  small  maximum  bank  angle.  The 
initial  point,  waypoints,  no-fly  zone,  and  target  from  Table  1  are  shown  in  Figure  13. 
This  figure  also  shows  the  baseline  trajectory  for  both  the  constant  speed  case  and 
decelerating  flight.  In  either  baseline  case,  the  control  law  is  to  make  a  maximum 
bank  turn  toward  the  next  waypoint  or  no-fly  zone  contact,  level  out,  and  then  fly 
straight  to  the  next  destination.  This  process  continues  until  target  intercept.  Later 
comparisons  will  illustrate  the  inefficient  overshoots  mentioned  earlier. 


Figure  13:  Baseline  Trajectories:  Constant  Speed  and  Decelerating  Flight 

4-2  2-D  Geometric 

A  geometric  approach  is  taken  for  the  special  case  of  the  HCV  at  constant  speed. 
A  geometric  technique  is  applicable  since  the  vehicle  traverses  a  constant  turn  radius 
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while  maintaining  a  constant  bank  angle.  A  bang-level-bang  controller  [23;  26;  88]  is 
optimal,  but  when  to  turn  remains  to  be  optimized.  Since  a  bang-level-bang  trajectory 
consist  only  of  straight  legs  and  turns,  waypoint  passage  must  occur  either  during  level 
flight  or  somewhere  in  a  turn.  If  the  waypoint  passage  occurs  along  a  straightaway, 
there’s  no  time  of  passage  to  optimize;  therefore,  the  point  of  waypoint  passage  along 
a  turn  is  what  is  to  be  optimized. 

Figure  14  shows  the  simplest  waypoint  intercept  case  involving  three  non-collinear 
points;  two  endpoints,  Pi  and  P2,  and  one  intermediate  waypoint,  Wp.  This  figure 
depicts  several  pilot  techniques  for  flying  from  P\ ,  directly  over  the  waypoint  Wp,  and 
then  onto  P2.  To  fly  the  profile  in  Figure  14(a),  the  pilot  aims  directly  for  the  way- 
point  and  upon  waypoint  passage  initiates  a  maximum  bank  turn  towards  endpoint 
P2.  This  is  the  baseline  technique  from  Section  4.1.  For  Figures  14(b)  and  14(c), 
the  pilot  initiates  a  heading  to  the  left  of  Wp,  then  upon  approaching  Wp  rolls  into 
a  maximum  bank  turn  which  is  timed  to  ensure  passage  over  the  waypoint  at  some 
fraction  along  the  turn,  and  finally  rolls  out  once  aligned  with  endpoint  P2.  The  last 
case  in  Figure  14(d),  where  the  pilot  initiates  a  heading  even  further  to  the  left  of 
Wp,  commences  a  maximum  bank  turn  such  that  waypoint  passage  occurs  exactly  at 
the  roll  out  to  endpoint  P2.  From  these  examples,  the  waypoint  passage  can  occur  at 
the  very  beginning  of  the  turn,  at  some  fraction  along  the  turn,  or  at  the  very  end 
of  the  turn.  The  fraction  of  turn  achieved  before  waypoint  passage  is  dictated  by  the 
orientation  of  the  turn  radius  with  respect  to  the  waypoint,  as  controlled  by  angle  y 
in  Figure  14;  therefore,  the  geometric  solution  aims  to  optimize  angle  y. 

For  this  analysis  a  constant  speed  is  assumed;  therefore,  minimizing  the  flight 
path  will  minimize  the  flight  time.  The  cost  becomes  the  flight  path  length  given  by: 

J  —  di{x)  +  RA9(x)  +  d2(x)  (44) 

Here  R  is  the  specified  turn  radius  and  the  values  for  d\,  d2,  and  A 9  are  all  geomet¬ 
rically  a  function  of  angle  y.  The  turn  radius  is  rotated  about  the  waypoint,  thus 
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(c)  x  =  0,  Optimal  Orientation  (d)  y  <  0,  Waypoint  on  Outgoing  Tangent 


Figure  14:  Turn  Radius  Orientation  Rotated  about  Intermediate  Waypoint  (Wp)  by 
Angle  X-  Radius  (R)  is  constant.  Endpoints  (Pi  and  P2),  and  Intermediate  Waypoint 
(Wp)  are  fixed. 

ensuring  passage.  Also,  the  initial  orientation  (x  =  0)  is  such  that  the  intersection  of 
the  tangent  lines  from  Pi  and  P2  is  collinear  with  the  waypoint  and  the  center  of  the 
turn  circle  as  shown  in  Figure  14(c).  The  angle  x  is  anchored  at  the  waypoint  and  is 
measured  counter-clockwise  to  the  center  of  the  turn  circle  as  shown  in  Figure  14.  A 
minimum  of  the  cost  J  with  respect  to  x  may  exist  when: 

dJ/dx  =  0  (45) 

Using  purely  geometric  relations,  the  solution  to  Equation  (45)  can  be  shown  to 
exist  at  y  =  0,  for  any  values  of  the  initial  setup  for  Pi,  P2,  Wp,  and  R.  This 
orientation  implies  that  waypoint  passage  occurs  at  the  halfway  point  along  the  turn, 
i.e.  the  midpoint  of  the  turn.  The  geometric  results  are  shown  in  Figure  15.  As 
expected,  graphically  the  optimal  geometric  solution  does  appear  shorter  than  the 
baseline  trajectory;  the  supporting  numerical  comparison  is  presented  in  Section  4.4. 
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60 1-  Geometric  (Constant  Speed) 


Baseline  (Constant  Speed) 
Baseline  (Decelerating) 


-40  -20  0 

Longitude  (deg) 


Figure  15:  Baseline  and  Geometric  Trajectories:  Constant  Speed  and  Decelerating 


The  significance  of  this  result  is  that  the  minimum  time  trajectory  can  be  ana¬ 
lytically  computed  for  any  given  initial,  final,  intermediate  waypoint,  and  turn  radius. 
Furthermore,  at  a  constant  velocity,  the  equations  of  motion  do  not  need  to  be  inte¬ 
grated,  so  long  as  a  relationship  between  maximum  bank  angle  and  minimum  turn 
radius  exist.  For  an  air  vehicle  to  produce  enough  lift  in  a  banked  turn  to  maintain 
level  flight,  the  centripetal  acceleration  ac  is  Lift  sin  a /m  where  Lift  must  be  mg/  cos  a. 
Thus  for  a  non-skidding  turn: 


V 2 

(46) 

ac  =  g  tan  cr  = 

~R 

V2 

(47) 

thus  R  =  - = 

Minimum  Turn  Radius 

g  tan  a 

In  this  geometric  approach,  the  no-fly  zones  can  be  treated  just  like  a  minimum  turn 
radius.  If  the  no-fly  zone  radius  is  greater  than  the  minimum  turn  radius,  a  tangent 
to  the  no-fly  zone  is  flown.  If  the  no-fly  zone  is  smaller  than  the  minimum  turn  radius, 
the  no-fly  zone  is  only  contacted  in  one  place  since  flight  along  the  boundary  cannot 
be  maintained. 

To  support  the  claim  of  midpoint  optimality,  the  optimal  configuration  shown 
in  Figure  14(c)  can  also  be  viewed  geometrically  as  a  circle  inscribed  within  a  vertex 
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of  a  triangle.  The  center  of  an  inscribed  circle  bisects  the  angle  of  each  vertex  [30]; 
therefore,  the  line  from  the  vertex  to  the  center  of  the  circle  is  the  angle  bisector. 
This  also  implies,  with  respect  to  the  center  of  the  circle,  that  the  angle  from  one 
tangent  to  the  angle  bisector  is  equal  to  the  angle  from  the  other  tangent  to  the  angle 
bisector.  Thus,  since  the  optimal  waypoint  in  Figure  14(c)  lies  on  the  angle  bisector, 
the  waypoint  occurs  halfway  along  the  arc  from  the  incoming  tangent  to  the  outgoing 
tangent.  The  optimality  of  this  halfway  point  is  proven  in  Appendix  C. 

For  an  application  with  multiple  waypoints  and  no-fly  zones,  a  simple  iterative 
technique  is  used  to  compute  the  optimal  geometry.  Initially  select  a  value  for  y, 
then  compute  the  point  of  intersection  of  the  two  tangent  lines  as  outlined  in  Ap¬ 
pendix  D.  This  provides  a  measure  as  to  whether  the  intersection  point  is  collinear 
with  the  waypoint  and  center  of  the  turn  circle.  This  is  an  iterative  process  since  a 
tangent  is  made  from  the  previous  and  following  turn  radii  or  no-fly  zones.  Since  the 
intermediate  waypoint  turn  radius  orientation  will  change  the  points  of  tangency,  the 
process  continues  until  no  more  changes  are  required.  Since  no  integration  is  required 
for  the  constant  speed  case,  each  iteration  is  fast,  and  convergence  requires  only  a  few 
iterations. 

4-3  2-D  Analytical,  Bryson’s  Method 

The  next  solution  method  is  an  analytical  approach  using  the  dynamic  opti¬ 
mization  problem  defined  in  Chapter  III.  This  section  references  Bryson’s  method 
because  in  Section  4.5  it  will  be  seen  that  the  purely  numerical  results  produce  dif¬ 
ferent  costate  histories  than  those  derived  here.  Studying  Bryson’s  method  provides 
insight  into  why  the  solution  must  hold  certain  properties,  this  intuition  is  lost  when 
the  numerical  method  in  Section  4.5  simply  produces  a  solution.  The  Breakwell  prob¬ 
lem,  presented  by  Bryson  [14:pg.l20]  and  Ross  et  al.  [78]  illustrates  the  differences 
on  a  simpler  problem  than  the  research  herein.  Reference  [78]  mentions  the  alternate 
form  of  the  Lagrangian  [40],  but  does  not  provide  an  explanation.  Section  4.6  explains 
the  fundamental  difference  in  problem  statement  and  compares  the  different  methods. 
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Appendix  F  specifically  addresses  the  Breakwell  problem.  It  will  be  shown  in  Sec¬ 
tion  4.6  that  Bryson’s  method  and  the  alternate  Covector/Costate  method  are  both 
correct;  however,  the  differences  must  be  understood  in  order  to  verify  the  numerical 
results.  Bryson’s  method  is  presented  first  because  it  provides  more  insight  into  the 
expected  solution. 

The  previous  chapter  provides  the  general  framework  for  the  dynamic  opti¬ 
mization  solution  for  the  problem  presented  herein.  In  order  to  solve  analytically, 
additional  steps  are  required.  The  vectors  u,  /x,  7rn,  and  7rm  are  required  in  order  to 
determine  the  optimal  control  and  to  propagate  the  costates  using  Equation  (12). 

4-3.1  Control.  The  Erst  task  is  to  determine  an  optimal  control  u.  The 
Hamiltonian  in  Equation  (9)  can  be  multiplied  out  to  show  that  there  is  only  one 
term  that  is  a  function  of  u  when  not  on  a  no-fly  zone,  i.e.  when  /x  =  0. 

H  =  \XV  cos  6  +  A yV  sin  9  +  A g - ^maxu  _|_  \ya  (43) 

Pontryagin  minimum  principle  yields  that  u  must  minimize  H  [2],  Furthermore,  since 
the  velocity  V  and  the  normalizing  parameter  tan  amax  are  always  positive,  u  must 
be  the  maximum  absolute  value  and  opposite  sign  of  A g  to  minimize  H .  This  defines 
A e  as  the  switching  function  for  u  in  a  bang-level-bang  controller  [18].  The  optimal 
control  specified  below  is  only  when  the  vehicle  is  not  on  the  no-fly  zone  constraint. 


(49) 


To  determine  the  control  while  on  a  no-fly  zone  constraint,  first  assume  that  the  no-fly 
zone  is  larger  than  the  minimum  turn  radius.  This  is  not  a  limitation,  it  just  implies 


This  u*  is  only  applicable  when 
off  the  no-fly  zone  constraints 


u*(t)  =  <  = 


1,  Xg  <  0 
0,  Xg  =  0 

-1,  Xg  >  0 
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that  if  the  no-fly  zone  is  smaller  than  the  minimum  turn  radius,  then  the  optimal 
control  would  still  be  Equation  (49).  The  tangent  entry  requirement  is  specified  as 
constraint  M  in  Equation  (28),  i.e.  S  —  0  and  S ^  =  0.  To  remain  on  the  no-fly 
zone  constraint,  the  S^2)  =  0  condition  must  remain  satisfied.  For  convenience,  Equa¬ 
tion  (26)  is  repeated  here  as  Equation  (50).  Starting  with  Equation  (50),  the  chain 
rule  and  state  derivatives  in  Equation  (22)  are  used  to  compute  the  time  derivatives 
of  S.  Setting  Equation  (52)  equal  to  zero,  the  optimal  control,  u*,  along  a  no-fly  zone 
boundary  is  computed  in  Equation  (53). 


5  = 

2  ( R j  ~  Axj  ~  Ay) ) 

(50) 

S W  = 

—A XjV  cos  9  —  AyjV  sin  9 

(51) 

s(2)  = 

— V2  +  u(Axj  sin  9  —  A yj  cos  6)  tan  amax 

+  (—A  Xj  cos  6  —  A  yj  sin  9)a  =  0 

(52) 

u*{t)  = 

V2 

This  u*  is  only  applicable  while 

(53) 

(A Xj  sin  9  —  Ai/j  cos  9)  tan  amax 

on  a  no-fly  zone  constraint 

The  term  in  front  of  a  in  Equation  (52)  is  zero,  since  it  is  a  multiple  of  Equation  (51), 
which  is  identically  zero  when  tangent  to  the  no-fly  zone  boundary.  It  is  not  eliminated 
from  Equation  (52)  since  the  partial  derivative  of  with  respect  to  the  state  vector 
is  still  required  to  propagate  the  costates. 

4-3.2  Costate  Propagation.  The  derivative  of  the  costates  must  be  deter¬ 
mined  in  order  to  propagate  the  costates  backwards  in  time.  The  costate  time  history 
is  used  to  compute  the  optimal  control  time  history.  Substituting  the  partial  deriva- 
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tive  of  Equations  (23)  and  (52)  into  the  transpose  of  Equation  (12)  produces: 


A  =  [-Apx]T-/,[S<2)]T 


Ax 

0 

A  y 

0 

Ae 

\XV  sin  9  —  XyV  cos  9 

•V 

\  a  \  ■  a  i  \  ^an 

—Ax  cos  9  —  Ay  sm  9  +  A g  — — ; — u 

L  y  Vz  -1 

u  sin  9  tan  amax  —  a  cos  9 
—u  cos  9  tan  crmax  —  a  sin  9 
(A Xj  sin  9  —  Ai/j  cos  9)  a 
-2V 


(54) 


(55) 


If  not  on  a  no-fly  zone  boundary,  then  /u  =  0  and  only  the  first  term  in  Equation  (55) 
needs  to  be  propagated.  However,  if  on  a  no-fly  zone  constraint,  then  /i  ^  0  and 
its  value  must  be  computed.  Solving  for  Hu  =  0  as  the  partial  derivative  of  H  in 
Equation  (19)  with  respect  to  u  yields: 


_  —  Ae 

^  V (A Xj  sin  9  —  Ayj  cos  9) 


(56) 


4-3.3  Final  Costates.  The  costates  are  initialized  at  the  final  time,  tf ,  to 
begin  the  backwards  in  time  integration.  The  vector  v  contains  two  elements,  vx 
and  uy,  since  there  are  two  elements  of  the  final  constraint  "0.  From  the  final  costate 
equation,  Equation  (13),  with  0X  =  0  and  0  only  a  function  of  x  and  y,  the  final 
costates  can  simply  be  defined  as: 


K(tf)  \{tf)  A e(tf )  A  v(tf) 


VX  Vy  ()  0 


(57) 


If  it  is  assumed  that  at  target  intercept  the  vehicle  is  not  in  a  turn,  then  u(tf)  =  0, 
A e(tf)  =  0,  and  A o(tf)  =  0.  Substitution  of  ux  and  uy  into  the  transversality  condition, 
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Equation  (15),  and  A e(tf)  =  0  into  Equation  (55),  leads  to: 


A  x{tf) 


—  cos  Of 


(58a) 


A y(tf) 


V,, 


y 


(58b) 


4-3-4  Jump  Conditions.  Next,  the  boundary  conditions  at  the  end  of  each 
segment  must  be  satisfied.  To  satisfy  the  intermediate  boundary  conditions  N  and 
M,  there  is  a  jump  in  the  Hamiltonian  and  the  costates  at  times  tt  and  tr  The 
jumps  occur  at  each  waypoint  passage,  and  entry  onto  the  tangent  of  each  no-fly  zone 
constraint. 

The  continuity  of  the  Hamiltonian  is  determined  by  the  jump  criteria  in  Equa¬ 
tion  (17b),  the  interior-point  constraint  N  for  the  waypoints,  and  the  interior-point 
constraint  M  for  the  no-fly  zones.  The  following  is  derived  from  the  definition  of  N 
in  Equation  (27)  and  from  the  definition  of  M  in  Equation  (28): 


ON/dU  =  0 
dM/dtj  =  0 


(59a) 

(59b) 


Since  the  jump  in  the  Hamiltonian  is  zero,  the  Hamiltonian  is  continuous  across  the 
interior-point  constraints: 


H(t~)  =  H(tt)  =  H(U) 
H(tj )  =  H(t+)  =  H(tj ) 


(60a) 

(60b) 


Also,  since  the  Hamiltonian  is  not  explicitly  a  function  of  time  it  must  be  a  constant; 
therefore,  solving  the  transversality  condition  in  Equation  (15)  and  substituting  A (tf) 
from  Equation  (13)  into  H  at  tf  in  Equation  (9),  the  constant  for  H  is  completely 
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determined  as: 


(61) 


H(tf)  =  H(t)  =  H(U)  =  H(tj)  =  -1 

The  costate  discontinuities  at  each  waypoint  are  determined  from  N  and  the 
jump  criteria  in  Equation  (17a).  The  control  is  continuous  if  waypoint  passage  at  tt 
occurs  along  the  turn,  i.e.  u{t~)  =  u(tf).  Since  the  assumption  is  that  the  waypoints 
are  sufficiently  spaced,  the  control  is  zero  just  before  starting  the  turn  at  time  tsi, 
u(t~i )  =  0,  and  non-zero  immediately  upon  commencing  the  turn,  zt(iT)  ^  0.  Similarly 
for  the  no-fly  zone,  the  control  is  zero  immediately  before  entering  the  no-fly  zone 
boundary,  u(tj)  =  0,  and  non-zero  immediately  after  contact,  u(tj~)  ^  0.  From  the 
waypoint  optimal  control,  u*(t),  defined  in  Equation  (49),  any  time  the  control  is  zero 
the  switching  costate  (A#)  is  zero.  As  derived  from  Equation  (17),  the  costate  jump 
at  each  waypoint  is  given  in  terms  of  the  unknown  constants,  irnx  and  i\ny\ 

WiV  =  M^+)T  +  [  Knx  nny  0  0]  (62) 

The  costate  propagation  in  Equation  (55)  shows  Xx  =  Xy  =  0;  therefore,  the  costates 
Xx  and  Xy  do  not  change  from  the  start  of  the  waypoint  turn  at  time  tsi,  to  immediately 
before  waypoint  passage  (referencing  forward  in  time): 

A  x(tsi)  =  Xx(t~)  =  Xx{tf)  +  nnx  (63a) 

A  y(tSi)  =  Xy(ti)  =  Xy(t+)  +nny  (63b) 

Using  Equation  (63)  and  the  continuity  of  the  Hamiltonian,  there  exists  a  relationship 
between  heading  at  the  start  of  the  turn,  0si,  and  the  heading  at  the  jump,  6t.  Using 
Equations  (55)  and  (63)  with  A  e(tSi)  =  0  and  A  (tf)  obtained  by  integrating  backwards 
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in  time  from  tf ,  the  derived  constants  are: 


7T  ri.T. 


-  sin  6i  (  A x(tf )  sin  9si  -  A y(tf )  cos  6si  ) 
cos (6i  -  6 si) 

cos  9 siXe{tl)  tan  amax  (u(tf)  -  u(t j-)) 


+ 


ft ny 


V"2 

v  si 

cos  Oi  (  A x(tf )  sin  9si  -  A y(tf )  cos  9si  ) 
cos(6f  -  9si) 

sin  9si\e(t+ )  tan  crmaa;  (u(tf )  -  u(t r )) 

V3 


(64a) 


(64b) 


When  waypoint  passage  occurs  along  the  turn  w(i^)  =  u(t)1")  and  only  the  hrst  terms 
in  Equation  (64)  are  non- zero.  Note  that  the  baseline  trajectory  defined  in  Section  4.1 
is  a  special  case,  where  the  passage  of  each  waypoint  occurs  at  the  start  of  the  turn, 
i.e.  tSi  =  ti  and  9si  =  9i,  but  the  turn  just  commenced  so  u(t~)  ^  u(fjt_);  thus,  the 
terms  at  the  end  of  Equation  (64)  are  non- zero.  Recall  the  baseline  is  a  simple  non- 
optimized  steer  and  point  control,  so  dynamic  optimization  is  not  required;  therefore, 
this  special  case  is  only  included  for  completeness  to  allow  the  Hamiltonian  to  be 
computed  for  the  baseline  to  see  how  it  compares  to  the  optimal  Hamiltonian. 

For  the  no-fly  zone  constraint,  the  components  of  M  in  Equation  (28)  must  be 
satisfied: 


s 

1  (jq  -  Ax]  -  Ay]) 

0 

sw 

—A XjV  cos  6  —  A yjV  sin  9 

0 

(65) 


The  control  is  zero  approaching  the  no-fly  zone,  which  implies  \g  has  been  zero, 
thus  the  entry  has  the  property  that  A $(tj)  =  0.  Also,  like  the  waypoint  jump, 
the  Hamiltonian  is  constant  H(tj)  =  H(t +).  Taking  the  partial  derivatives  of  M 
in  Equation  (28)  and  solving  for  the  constants  using  A o(tj)  =  0,  A o(tj)  =  0,  and 
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Hit- )  =  Hit t)  leads  to  the  following: 


Axj  =  x(tj )  -  ,  A yj  =  y(tj )  -  ycj 


K(tj ) 

K(tj  )  -  KmoAxj  -  KmlVj  cos  9j 

) 

\y(t+ )  -  TTmoAyj  -  'KmlVj  sin  9j 

A  e(tj) 

A e(t^ )  +  7rmi  (A Xj  sin  9j  -  Ayj  cos  9j)  Vj 

_  MtJ)  _ 

_  Ay(t+) 

A^fb )  sin  —  Ay(td" )  cos  9j 

Axj  sin  9j  —  Ai/j  cos  9j 
\g(tf)u(tf)  tan  amax  _  -Afl(t+) 

1A3  Vj  (A Xj  sin  9j  —  Ay3  cos  9j) 


(66) 


(67) 


(68) 

(69) 


Since  the  states  at  time  tj  are  continuous,  the  plus  and  minus  superscripts  are  unnec¬ 
essary.  The  control  u(tij~)  in  Equation  (69)  is  replaced  with  the  u  required  to  maintain 
=  0  from  Equation  (53).  Within  the  derivation  of  Equation  (69),  additional  terms 
appear  that  are  multiples  of  S'd);  however,  since  these  terms  are  not  shown 

for  simplification. 

Having  now  defined  all  the  criteria  for  the  optimal  trajectory,  the  next  step  is 
to  connect  all  the  solution  pieces.  The  geometric  solution  provides  an  initial  guess  of 
the  ten  solution  times;  the  waypoint  passage  times  (2),  the  no-fly  zone  entry  time  (1), 
the  final  time  (1),  the  end  of  the  initial  turn  toward  the  first  waypoint  (1),  the  start 
and  end  time  of  each  waypoint  turn  (4),  and  the  no-fly  zone  exit  time  (1).  These 
times,  and  the  associated  optimal  control,  are  used  to  propagate  the  states  forward. 
From  the  forward  integration,  the  accuracy  of  hitting  the  final  target,  Equation  (25); 
passing  through  each  waypoint,  Equation  (27);  and  satisfying  the  tangency  conditions, 
Equation  (28),  are  all  used  as  the  measure  of  solution  success.  At  this  point,  for  the 
mission  as  defined  in  Table  1,  the  solution  has  eight  solution  criteria.  The  forward 
states  are  interpolated  during  the  backward  integration  for  the  costates.  During  the 
backward  integration,  the  appropriate  derivatives  and  jumps  are  computed.  The 
accuracy  of  achieving  each  A e(tjj)  =  0  are  additional  solution  criteria;  two  for  this 
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mission  since  there  are  two  waypoints.  This  increases  the  total  number  of  solution 
criteria  to  ten.  Thus,  a  search  vector  of  ten  times  is  used  to  satisfy  ten  success  criteria, 
which  should  each  be  zero.  A  root  solver  is  used  to  manipulate  a  search  vector  to  drive 
the  solution  vector  to  zeros.  For  this  research,  the  solution  or  search  times  are  input 
into  the  Matlab®  root  solver,  f  solve,  to  drive  the  solution  criteria  towards  zero.  A 
related  indication  of  success  (optimality)  is  to  observe  the  value  of  the  Hamiltonian. 
During  this  derivation,  H(t)  is  assumed  continuous,  but  only  the  optimal  solution 
will  yield  H(t)  =  —1  from  Equation  (61)  for  all  t.  The  results  from  this  process  are 
presented  in  the  next  section. 

4-4  2-D  Comparison,  Geometric  versus  Analytical 

This  comparison  includes  the  baseline,  geometric,  and  analytical  dynamic  op¬ 
timization  approach  for  constant  speed  and  deceleration.  It  is  presented  prior  to 
the  numerical  results  because  the  numerical  results  converge  to  the  same  dynamic 
optimization  solution,  just  the  technique  has  changed.  This  section  will  therefore 
illustrate  the  time  savings  possible  using  optimal  control.  Several  HCV  trajectories, 
based  on  the  mission  described  in  Table  1,  are  presented  in  Figure  16. 


Figure  16:  Baseline,  Geometric,  and  Analytical  Dynamic  Optimization  Trajecto¬ 

ries;  both  Constant  Speed  and  Decelerating 
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For  this  problem  the  initial  heading  is  fixed  and  the  final  heading  is  not  con¬ 
strained;  therefore,  the  optimization  is  the  orientation  angle  for  the  turn  radii  of 
the  first  and  second  intermediate  waypoints,  and  the  point  entering  the  no-fly  zone 
boundary.  The  baseline,  geometric,  and  dynamic  optimization  numerical  results  are 
presented  in  Table  3.  For  constant  speed,  the  geometric  and  dynamic  optimization 
results  are  identical,  as  derived  in  Appendix  C.  With  deceleration,  the  geometric  and 
dynamic  optimization  results  are  nearly  the  same,  and  are  thus  indistinguishable  in 
Figure  16. 

The  control  u  is  bang- level-bang  based  on  the  value  of  the  costate  A g.  The 
costates  and  control  are  plotted  in  Figure  17(a).  Throughout  the  rest  of  this  docu¬ 
ment,  the  labeling  and  scaling  for  some  variables  may  be  along  the  right  side  of  the 
figure,  e.g.  Xx,  A g,  and  u  in  Figure  17(a).  In  Figure  17(a)  the  value  of  u  switches  with 
the  sign  of  A g.  The  exception  to  the  bang- level- bang  control  is  while  the  trajectory 
is  on  the  no-fly  zone  boundary.  The  value  of  u  is  less  than  the  maximum  to  avoid 
violating  the  constraint;  however,  optimality  is  maintained  since  Hu  =  0  is  maintained 
by  Equation  (56). 

In  Figure  17(b)  the  effect  of  a  non-optimal  solution  is  demonstrated  by  observ¬ 
ing  the  value  of  the  Hamiltonian.  Here,  the  optimal  solution  correctly  maintains  a 
Hamiltonian  value  of  —1;  however,  with  deceleration  the  geometric  solution  is  seen 
to  be  non-optimal.  Non-optimality  is  expected  with  the  geometric  solution  since  it 
forces  constant  turn  radii;  however,  this  no  longer  represents  a  maximum  bank  angle 
turn  while  decelerating.  The  baseline  Hamiltonian,  not  plotted,  is  even  further  from 
optimal  reaching  a  Hamiltonian  value  of  1.04  (that  is  2.04  from  the  optimal  value  of 
-!)• 
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(a)  Costates  and  Control  u  Time  History  for  Dynamic  Optimal  Control  of  Decelerating  Flight 
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(b)  Optimal  Hamiltonian  H(t)  =  —1,  and  Non-Optimal  H(t )  ^  —1. 

Figure  17:  Nondimensional  Costates,  Control,  and  Hamiltonian  Time  Histories 


In  Table  3  there  are  some  notable  trends  beyond  just  minimizing  flight  time.  For 
the  constant  speed  baseline  case,  the  vehicle  does  not  meet  the  design  specifications 
of  9,000  nmi  in  2  hrs.  So,  not  only  does  the  baseline  take  14%  longer,  it  fails  the 
mission.  A  mission  failure  may  be  characterized  by  consuming  all  fuel  prior  to  reaching 
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Table  3:  Numerical  Results  for  Each  Trajectory  9,000  nmi  in  2  hrs  at  100,000  ft 


Vmin  =  310  KEAS,  Initial  Conditions:  do  =  0  deg,  M0  =  7.66 
V0  =  2.315  km/s,  r0  =  6408.48  km,  g0  =  9.7169  m/s2 
Technique  Decel  (m/s2)  Vf  KEAS  Mf  Distance  (nmi)  Time  (hlumnuss) 


Baseline 

0 

531.69 

7.66 

10285 

02:17:08 

Geometric 

0 

531.69 

7.66 

9000 

02:00:00 

Optimal 

0 

531.69 

7.66 

9000 

02:00:00 

Baseline1 

0.1 

301.31 

4.34 

9822 

02:47:11 

Geometric 

0.1 

331.24 

4.78 

8853 

02:25:28 

Optimal 

0.1 

331.89 

4.78 

8831 

02:24:59 

1  Invalid  solution 

since  minimum  airspeed  for  level  flight  ( Knm 

)  is  not  maintained. 

Table  4: 

Nondimensioual  Times  for  Turns  and  Target  Arrival 

Initial 

Waypoint  1 

Waypoint  2 

No-Fly  Zone 

Target 

Technique 

end 

start 

pass  end 

start 

pass  end 

enter  exit 

arrive 

Baseline 

0.507 

2.616 

2.616  4.339 

5.819 

5.819  7.228 

8.328  9.157 

10.132 

Geometric 

0.688 

1.843 

2.738  3.633 

4.487 

5.415  6.344 

6.984  7.894 

8.866 

Optimal 

0.688 

1.843 

2.738  3.633 

4.487 

5.415  6.344 

6.984  7.894 

8.866 

Baseline 

0.502 

2.748 

2.748  4.149 

6.323 

6.323  7.099 

9.031  10.713 

12.352 

Geometric 

0.641 

2.084 

2.842  3.624 

5.276 

5.940  6.623 

8.129  9.250 

10.747 

Optimal 

0.615 

2.114 

2.834  3.542 

5.288 

5.918  6.539 

8.116  9.218 

10.712 

All  trajectories  to  = 

0.  First 

;  three  entries 

are  constant  speed,  last  three  are  decelerating. 

the  target.  A  similar  time  savings  is  seen  with  the  decelerating  case;  however,  the 
baseline  is  unsatisfactory  since  it  has  reached  an  airspeed  too  slow  to  sustain  level 
flight  (Vmin  =  310  KEAS  for  level  flight),  thus  this  trajectory  is  invalid.  Once  again, 
the  optimized  trajectory  has  not  only  saved  time,  but  enabled  mission  success.  Table 
4  is  provided  as  a  means  of  verifying  duplication  of  the  presented  results. 


4-5  2-D  Numerical 

This  section  is  an  introduction  to  the  numerical  technique  required  to  solve 
the  more  complex  3-D  CAV  problem.  These  numerical  results  are  computed  using 
the  software  package  GPOCS,  and  are  compared  to  the  analytical  results  previously 
computed  in  Section  4.3. 
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Figure  18  shows  that  the  numerical  results  are  approximately  equal  to  the  ana¬ 
lytically  derived  results.  It  also  illustrates  the  numerical  technique  of  separate  phases, 
where  the  waypoints  and  target  are  specified  as  terminal  boundary  conditions  for  the 
appropriate  phase. 


(a)  Map 


(b)  Control 

Figure  18:  Analytical  (Bryson’s  method)  versus  Numerical  Results 


The  costates  in  Figure  19  match  for  most  of  the  time  history;  however,  the  nu¬ 
merical  results  have  different  jump  values  and  derivatives  while  the  trajectory  is  on  the 
no-fly  constraint.  In  comparison  to  Bryson’s  method,  the  GPOCS  software  package 
uses  an  alternate  path  constraint  method  that  is  compatible  with  the  Gauss  pseu- 
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dospectral  costate  mapping  theorem  [4],  shown  in  Figure  11(b).  To  understand  the 
differences,  this  alternate  method  is  used  to  derive  another  set  of  analytical  equations 
for  the  costates;  the  derivation  and  comparison  is  presented  next  in  Section  4.6. 


Figure  19:  Costates:  Analytical  (Bryson’s  method)  versus  Numerical  Results 


4-6  2-D  Analytical,  Alternate  Method 

The  Gauss  pseudospectral  costate  mapping  theorem  is  used  within  GPOCS  to 
convert  the  discrete  KKT  multipliers  into  the  continuous  costates.  The  user  supplied 
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path  function  provided  to  GPOCS  is  simply  S,  the  no-fly  zone  constraint.  Tims  S, 
not  S(2\  is  adjoined  to  the  Hamiltonian: 

H  =  L  +  XTf  +  fj,S  (70) 

This  change  in  Hamiltonian  implies  a  change  in  the  costate  propagation  equations, 
as  compared  to  Equation  (54): 

A  =  [-A  Tfx]T-»[Sx]T 


^ X 

0 

1 

1 

> 

Xy 

0 

= 

-  h 

Xe 

XxV  sin  9  —  XyV  cos  9 

0 

Ay 

\  A  \  ■  A  1  \ 

— Xx  cos  9  —  Xy  sm  9  +  Xg — — — u 

0 

This  change  in  adjoined  path  constraint  has  created  non-zero  derivatives  for  Xx  and  \y. 
Furthermore,  the  previously  derived  interior-point  constraint,  M,  is  not  imposed  so 
there  is  no  jump  in  the  costates.  From  the  new  definition  of  A g  in  Equation  (71)  there 
is  zero  change  in  its  derivative  occurring  at  no-fly  zone  contact;  therefore,  without  a 
change  in  derivative  it  must  remain  its  previous  value  of  Ag  =  0.  This  criteria  leads 
to  the  following  relationship  along  the  no-fly  zone  boundary: 

\XV  sin  9  —  XyV  cos  9  —  0  (72) 

Since  V  is  always  positive  it  can  be  divided  out  of  this  equation.  The  next  step  is  to 
take  the  time  derivative  of  both  sides  of  the  equation: 

sin  6Xx  +  A*  cos  9  9  —  cos  9Xy  +  Xy  sin  9  9  =  0  (73) 
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Substituting  in  the  values  of  =  fiAxj  and  Xy  =  yAyj,  from  Equation  (71),  leads 
to  an  expression  for  //: 

_  ~(Xx  cos  9  +  Xy  sin  9)  V 
^  (A Xj  sin  9  —  Ay3  cos  0)2 

Upon  no-fly  zone  contact  the  costates  are  propagated  using  Equations  (71)  and  (74)  to 
create  the  time  history  presented  in  Figure  20.  The  control  time  history  is  unchanged 
since  u*(t)  in  Equation  (53)  is  not  a  function  of  the  costates. 


Figure  20:  Analytical  (Alternate  Method)  versus  Numerical  Results 
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Figure  21:  Numerical  and  Analytical  Path  Multiplier,  /i 

The  correlation  in  the  numerical  results,  as  shown  in  Figure  20,  with  the  al¬ 
ternate  derivation  of  the  costate  propagation  equations  clearly  shows  that  numerical 
computation  is  conforming  to  the  optimality  criteria.  This  understanding  of  what 
path  constraint  is  being  adjoined  to  the  Hamiltonian  to  create  the  time  history  of  the 
costates  is  used  in  the  CAV  3-D  analysis.  Part  of  the  numerical  solution  by  GPOCS 
is  the  value  of  the  path  multiplier,  /i.  The  value  of  y  from  Equation  (74)  is  plotted 
in  Figure  21,  along  with  the  GPOCS  path  constraint  multiplier  output.  The  GPOCS 
values  have  not  been  altered;  however,  the  computed  y  values  are  set  to  zero  for  all 
times  off  the  no-fly  zone  boundary  constraint.  The  3-D  model  is  analyzed  next  using 
this  understanding  gained  from  the  simpler  2-D  results. 

4-7  3-D  Analytical 

The  analytical  solution  is  derived  first  to  determine  the  form  of  the  optimal 
control.  The  first  step  is  to  derive  the  solution  with  no  path  constraints.  Then,  with 
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path  constraints,  determine  how  the  propagation  of  the  costates  changes  with  the 
inclusion  of  the  path  multipliers,  /iT. 

4-7.1  3-D  Optimal  Unconstrained  Control.  First  assume  the  Hamiltonian  is 

only  a  function  of  the  integrand  and  the  costates,  i.e.  p  =  0  since  the  path  constraints 
are  assumed  inactive: 


H  =  L  +  XTf 

=  XxV  cos  6  +  XyV  sin  9  +  XhV'y  —  Ay 


+  X~, 


BVe  dr°hC£  cos  cr  —  —  +  V 


BV2e-Pr°h(l  +  c2£) 

2  E* 

+  XeBVe~0rohcesma 


(75) 


The  first-order  optimality  condition  requires  Hu  =  0  and  the  second-order  criteria 
requires  Huu  >  0  [12:pg.66;  50],  i.e.  positive  definite.  These  represent  the  following 
matrices: 


u 


Hn 


H 


UU 


T 


cr  ci 

Ha 

Hce 

= 

Haa 

Hace 

HCia 

HcfCf 

(76) 

(77) 

(78) 


Starting  with  Equation  (75)  the  partial  derivative  of  the  Hamiltonian  with  respect  to 
the  controls  are: 


Ha  =  -X1BVe-f}rohcesma  +  XgBVe-^ce  cos  a  =  0  (79a) 

nv2p-Pr°h 

HCl  =  -Xv - - - -\-  XryBV e~^r°h  cos  u  T  XgBV e~^r°h  sin  u  =  0  (79b) 
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The  second  partial  derivative  is  taken  to  ensure  Huu  is  positive  definite: 


H 


UU 


—ciBVe  /3r°/l(A7  coscr  +  coscx)  BVe  A7  sin  a  +  \g  cos  a) 

BV2e~/3r°h 

BVe~Pr°h(— A7  sin  a  +  A#  cos  a)  —Ay - — - 


>  0 
(80) 


The  off  diagonals  of  Huu  are  zero  since  they  are  simply  H^/ce,  where  Ha  =  0  from 
Equation  (79a).  Thus,  the  eigenvalues  of  Huu,  which  must  be  >  0  [12:pg.66],  are 
equal  to  the  terms  on  the  main  diagonal: 


—cgBVe  /3r°/l(A7  cos  a  +  Ag  cos  a)  >  0 
,  BV2e~f5r°h 


In  analyzing  Equation  (81b);  the  terms  B ,  V2,  e  @ r°h ,  and  E*  are  all  positive,  therefore 
Ay  must  be  negative: 

Ay  <  0  (82) 

In  order  to  satisfy  the  conditions  in  Equation  (79),  two  equations  must  be  solved 
simultaneously: 


(81a) 

(81b) 


— A7  sin  a  +  \q  cos  a  =  0  (83a) 

Ay cnjE*  +  A7  cos  a  +  \q  sin  a  =  0  (83b) 

To  solve  Equation  (83a),  assume  A7  =  kcosa  and  A g  =  k  sin  a  where  k  is  a  constant 

to  be  determined.  From  the  assumed  values  for  A7  and  A#,  and  the  trigonometric 

identity  sin2  a  +  cos2  a  =  1  [51:pg.A52],  the  following  must  also  be  true: 

A2  +  A2  =  k2  (84) 
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Thus,  an  intermediate  solution  is: 


±^7  +  ^ 

(85a) 

±cosa^X2  +  Ag 

implies 

cos  cr  =  ±A7/ y/ X2+X2g 

(85b) 

±  sin  a  ■J  X2  +  \2e 

implies 

sin  cr  =  ±\g/  J  X2  +  X2e 

(85c) 

The  sign  of  k  must  still  be  determined  from  Equation  (81a);  enforcing  the  control  lim¬ 
itation  Cf  >  0,  implies  A7  cos  cr  +  Ag  cos  cr  <  0.  Using  the  solution  from  Equation  (85) 
leads  to: 

(±A?  ±  AD/^  +  AI  <  0  (86) 

To  satisfy  Equation  (86)  implies  k  must  be  negative;  therefore,  the  solution  for  the 
control  a  is: 


cos  a  = 

-A7/ ^X2  +  X2e 

(87a) 

sin  a  = 

-Xg/  ^  X2  +  Ag 

(87b) 

a  = 

atan2(  sin  cr,  cos  cr) 

(87c) 

Equation  (87)  is  now  plugged  into  Equation  (83b)  to  solve  for  q  as: 


Ci 


VM 


(88) 


The  atan2  function  in  Equation  (87c)  solves  for  the  correct  quadrant  for  the  angle  cr. 
This  may  seem  unnecessary  since  a  is  confined  to  —60°  <  a  <  60°;  however,  the  sign 
of  a  is  still  determined  by  the  unbounded  solution.  Thus,  using  this  methodology 
maintains  the  correct  sign  for  cr,  and  avoids  incorrectly  jumping  from  one  bound  to 
the  next  when  the  unbounded  solution  exceeds  ±90°. 
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The  numerically  solved  values,  from  the  next  section,  for  Ay,  A7,  and  A g  will 
be  input  into  Equations  (87)  and  (88)  to  integrate  the  states  forward  to  verify  the 
numerical  results  are  adhering  to  the  optimality  criteria  enforced  in  this  derivation. 

4-7.2  3-D  Optimal  Constrained  Control.  The  CAV  model  allows  for  two 

controls,  a  and  q,  which  have  specified  limitations.  It  will  be  shown  that  when  one 
control  is  restricted  to  its  limit,  the  other  control  law  must  be  modified  to  maintain 
optimality.  In  other  words,  if  a  control  is  at  its  limit,  the  control  laws  derived  in 
Equations  (87)  and  (88)  may  no  longer  be  valid. 

The  control  limitations,  Equation  (33),  are  adjoined  to  the  Hamiltonian  via  the 
multiplier  /i,  as  defined  in  Equation  (9),  repeated  here: 

H  =  L  +  XTf  +  \iTC  (89) 

As  an  example,  assume  the  bank  angle  has  reached  its  maximum  positive  value,  thus 
the  bank  angle  constraint  is  active.  For  optimality,  Ha  and  HCe  must  still  equal  zero; 
however,  since  fi  is  now  non-zero  the  solutions  in  Equations  (87)  and  (88)  are  no 
longer  valid.  The  new  equations  for  Ha  and  HCf  are: 

Ha  =  -BVe-prohXJcea  +  BVe-l3rohXececosa  +  ia  =  0  (90a) 

BV2p-Pr°hr<, 

HCl  =  Ay - - - -  +  X1BV e~Pr°h  cos  a  +  XqBV e~l3r°h  sin  a  —  0  (90b) 

E* 

Unlike  the  previous  section,  a  is  now  known  to  be  at  its  maximum 
angle  is  at  its  maximum  positive  value  the  optimal  q  is: 

(A7  cos  a  +  Xq  sin  a)E* 


Thus,  when  bank 


(91) 
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x  10“3 


Figure  22:  Control  Multiplier  /i  for  Constrained  Bank  Angle  (er  =  60°) 


To  compare  this  derivation  to  the  numerical  output,  Equation  (91)  is  used  in  Equa¬ 
tion  (90a)  to  define  //: 


fi  =  BVe  /3r°h(A7  sin  a  —  \e  cos  <t)q  (92) 

The  comparison  between  the  numerically  generated  results  and  this  analytical  form 
is  presented  in  Figure  22.  The  time  span  in  Figure  22  is  at  the  end  of  the  trajectory 
where  maximum  bank  angle  is  required  to  achieve  the  optimal  solution.  The  matching 
results  demonstrate  that  the  numerical  results  are  satisfying  the  optimality  criteria. 
For  the  case  of  maximum  q,  the  optimality  criteria  in  Equations  (87)  and  (88)  remain 
valid  since  Equation  (79a)  is  not  influenced  by  q  (q  is  present  but  it  is  divided  out 
thus  leaving  the  originally  derived  solution). 
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Figure  23:  Map:  Seven  Phase  Numerical  Results 

4-8  3-D  Numerical  Results 

To  maximize  the  accuracy  of  the  numerical  results  the  trajectory  is  broken  up 
such  that  every  jump  in  the  costates  occurs  at  the  end  of  a  phase.  This  requires  an 
iterative  approach  to  produce  these  more  accurate  results;  however,  the  more  accurate 
results  are  used  to  verify  optimality  and  to  test  the  solution  convergence  with  fewer 
phases,  i.e.  less  iterations  from  previous  solutions. 

The  costate  jumps  occur  at  contact  with  the  heating  and  no-fly  zone  constraints. 
Preliminary  numerical  results,  with  fewer  phases,  are  used  to  identify  the  number  of 
path  constraint  contacts.  Next,  additional  phases  are  specified  to  more  accurately 
capture  these  events  and  verify  the  original  results.  Thus,  this  trajectory  is  broken  into 
seven  phases;  to  to  heating  constraint  contact,  heating  constraint  contact  to  waypoint 
1,  waypoint  1  to  no-fly  zone  1,  no-fly  zone  1  to  waypoint  2,  waypoint  2  to  first  contact 
of  no-fly  zone  2,  first  contact  of  no-fly  zone  2  to  second  contact  of  no-fly  zone  2, 
second  contact  of  no-fly  zone  2  to  target,  as  seen  on  the  map  in  Figure  23.  Figure  24 
shows  the  states  throughout  the  seven  phases.  One  important  factor  is  that  the 
states  remain  continuous  across  the  phase  breaks.  These  numerical  results  are  using 
collocation  techniques,  not  integration,  so  verifying  the  results  satisfy  the  equations 
of  motion  is  presented  in  the  next  section.  The  magnitude  of  flight  path  angle,  7, 
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V  (km  /  sec) 


Figure  24:  States:  Seven  Phase  Numerical  Results 
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7  (deg) 


remains  small  throughout  the  trajectory;  therefore,  the  small  angle  approximation  is 
still  valid.  Also,  the  specified  final  altitude  is  seen  to  be  satisfied  in  the  top  graph  of 
Figure  24. 

The  costates  are  presented  in  Figure  25.  The  jumps  in  the  costates  are  at  the 
end  of  each  phase,  as  expected.  The  zero  derivatives  of  Xx  and  Xy  are  clearly  seen.  The 
terminal  boundary  conditions  for  Ay,  A7,  and  Xq  are  all  met  by  equalling  zero,  thus 
satisfying  Equation  (13)  using  Equation  (34).  Finally,  the  Hamiltonian  is  constant 
and  continues  to  satisfy  H(t)  =  —1. 

The  controls  are  plotted  in  Figure  26.  The  controls  are  seen  to  remain  within 
their  specified  limitations;  for  example,  q  remains  on  the  limit  of  q  =  2  at  ap¬ 
proximately  200  seconds  into  the  trajectory,  and  bank  angle  remains  at  the  limit 
of  a  =  60°  at  approximately  2700  seconds.  The  path  constraint  multiplier  for  this 
maximum  bank  angle  constraint  contact  was  previously  plotted  in  Figure  22. 

The  path  constraints  are  shown  in  Figure  27,  with  the  no-fly  zone  constraints  S 
labeled  on  the  left  and  the  heating  constraint  Q  labeled  on  the  right.  Expanded  views 
of  these  path  constraints  are  plotted  in  Figure  28.  Contact  with  a  path  constraint 
occurs  when  the  value  of  the  constraint  equals  zero,  i.e.  S  =  0  from  Equation  (35)  or 
Q  =  0  from  Equation  (41).  For  each  path  constraint  there  is  contact  without  travel 
along  the  constrained  arc.  In  Figure  28(c)  there  are  two  separate  contacts  on  no-fly 
zone  2,  hence,  the  two  sets  of  costate  jumps  labeled  in  Figure  25. 

The  next  section  takes  the  numerical  costate  histories  in  Figure  25,  computes  the 
analytical  controls  from  Sections  4.7.1  and  4.7.2,  and  integrates  the  states.  The  intent 
is  to  verify  that  the  numerical  results  are  properly  estimating  the  optimal  solution. 
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Figure  25:  Costates:  Seven  Phase  Numerical  Results 
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S,  No-Fly  Zone  1  (km/1e7) 


Figure  26:  Controls:  Seven  Phase  Numerical  Results 


Figure  27:  Path  Constraints:  Seven  Phase  Numerical  Results 
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Figure  28:  Path  Constraints  Expanded:  Seven  Phase  Numerical  Results 


4-9  3-D  Comparison,  Analytical  versus  Numerical 

In  Figures  29  through  33  the  numerical  results  are  compared  to  the  form  of  the 
solution  analytically  derived  in  Sections  4.7.1  and  4.7.2.  The  states  and  constraints 
are  nearly  a  perfect  match.  For  the  controls  in  Figure  31  the  analytical  results  appear 
smoother,  especially  at  the  beginning  of  the  trajectory.  This  is  attributed  to  costate 
interpolation  used  to  compute  the  controls.  The  derivations  in  Section  4.7.2  discussed 
the  change  in  optimal  control  law  when  a  control  constraint  is  reached.  The  correlation 
in  the  results  at  the  end  of  the  trajectory  in  Figure  30  demonstrate  the  successful 
implementation  of  the  revised  control  law.  The  data  point  jump  at  the  last  point 
in  Figure  31  is  caused  from  the  singularities  that  arise  as  the  costates  Xy,  A7,  and 
A e  approach  the  analytically  derived  values  of  zero.  Lastly,  in  the  expanded  view  in 
Figure  33  there  are  differences  in  constraint  contact  points.  Herein  lies  the  instability 
with  the  integration  in  the  shooting  method,  the  points  of  contact  and  departure  must 
be  exact  or  the  propagated  costates  may  diverge. 
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Figure  29:  Map:  Seven  Phase  Numerical  and  Analytical 


Figure  30:  States:  Seven  Phase  Numerical  and  Analytical 
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Figure  31:  Controls:  Seven  Phase  Numerical  and  Analytical 
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Figure  32:  Path  Constraints:  Seven  Phase  Numerical  and  Analytical 
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Figure  33:  Path  Constraints  Expanded:  Seven  Phase  Numerical  and  Analytical 


4-10  Numerical  Comparison  of  Phase  Breakpoints 

This  rigorous  analysis  forces  this  problem  to  be  broken  up  into  many  phases, 
several  required  advanced  knowledge  of  the  solution,  i.e.  the  number  of  boundary 
contacts.  The  next  objective  is  to  evaluate  the  solution  effectiveness  using  fewer 
phases.  The  map  results  are  seen  in  Figure  34. 


Figure  34:  Map:  Number  of  Phases  Comparison 

For  this  mission,  there  are  two  known  intermediate  waypoints;  therefore,  there 
are  at  least  three  phases;  to  to  U= i,  ti= i  to  U= 2,  and  U= 2  to  tf.  Through  some  solution 
knowledge,  but  without  confining  the  solution,  a  fourth  phase  is  added.  This  phase 
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ends  at  the  first  point  the  flight  path  angle  passes  through  zero.  The  rational  is  the 
first  assault  on  the  heating  limit  requires  the  most  accuracy;  therefore,  by  placing  a 
phase  transition  here,  near  the  potential  heating  constraint  contact,  additional  solu¬ 
tion  accuracy  is  achieved.  The  rigorous  seven  phase  results  are  compared  to  the  much 
simpler  four  phase  results  in  Figures  34  through  39. 


Figure  35:  States:  Number  of  Phases  Comparison 
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Figure  36:  Costates:  Number  of  Phases  Comparison 
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Figure  37:  Controls:  Number  of  Phases  Comparison 


Figure  38:  Path  Constraints:  Number  of  Phases  Comparison 
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Figure  39:  Path  Constraints  Expanded:  Number  of  Phases  Comparison 


The  closeness  in  these  results  demonstrates  that  the  jump  in  costates  due  to  path 
constraint  contact  can  be  estimated  as  continuous.  Quantitatively,  the  final  times  for 
the  seven  phase  and  four  phase  solutions  are  2719.42  sec  and  2719.18  sec,  respectively; 
thus,  well  within  the  accuracy  of  the  model.  Future  analysis  is  required  to  generalize 
this  small  error  to  higher  numbers  of  waypoints,  no-fly  zones,  and  heating  constraint 
contacts.  The  setup  for  the  four  phase  problem  is  much  simpler  and  it  does  not  require 
advanced  knowledge  of  the  number  of  times  the  path  constraints  are  contacted.  Due 
to  the  direct  connection  between  phases  and  waypoint,  regardless  of  heating  or  no-fly 
zones,  it  makes  it  possible  to  automate  the  setup  based  on  the  mission  objectives. 


4-11  Summary  of  Analysis 

The  development  herein  built  upon  itself  to  eventually  determine  and  verify  a 
viable  optimal  reentry  trajectory  solution  technique.  The  simpler  2-D  HCV  case  is 
used  as  a  foundation  to  understand  and  introduce  the  analytical  and  numerical  tech¬ 
niques.  The  geometric  solution  provides  an  independently  developed  optimal  solution 
to  validate  both  the  analytical  and  numerical  solutions.  Bryson’s  path  constraint 
method  is  seen  to  differ  from  that  used  in  the  numerical  collocation  software.  Un¬ 
derstanding  this  difference  provides  insight  into  which  costates  are  expected  to  have 
discontinuities;  therefore,  provides  a  means  of  replicating  the  results  using  analytical 
methods.  The  more  complex  3-D  CAV  model  demonstrates  the  ability  to  incorporate 
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all  of  the  vehicle  and  mission  constraints.  Using  dynamic  optimal  control  theory  an 
analytical  solution  is  derived  to  validate  these  numerical  results.  After  verifying  that 
the  numerical  results  are  indeed  satisfying  the  optimality  criteria,  the  next  step  is  to 
simplify  implementation.  By  replicating  the  numerical  results  using  just  four  phases 
instead  of  the  full  seven  phases  demonstrates  that  a  priori  knowledge  of  path  con¬ 
straint  contact  is  unnecessary;  therefore,  an  automated  problem  setup  technique  can 
be  developed  primarily  dependent  on  the  number  of  waypoints.  The  following  section 
explores  the  importance  of  these  results. 
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V.  Conclusions,  Contributions,  and  Future  Work 

Having  presented  the  problem,  analysis,  and  results  the  next  step  is  to  address 
the  conclusions  to  be  drawn,  how  they  affect  the  user  and  the  research  community, 
and  what  future  work  lies  ahead. 

5. 1  Conclusions 

This  research  successfully  determined  and  verified  a  numerical  technique  capable 
of  generating  an  optimal  reentry  trajectory  satisfying  all  of  the  specified  vehicle  and 
mission  constraints.  The  suitability  for  satisfying  the  Global  Strike  mission  also  in¬ 
cludes  the  autonomy  and  user  implementation  of  the  solution  technique.  The  solution 
is  considered  autonomous  due  to  several  factors;  manual  manipulation  of  variables  is 
not  required,  the  problem  has  a  repeatable  setup  structure,  and  the  solution  can  con¬ 
verge  even  with  a  very  simple  “point-to-target”  initial  guess.  The  repeatability  is  the 
structured  setup  based  on  the  number  of  waypoints  and  convergence  independent  of 
the  number  of  no-fly  zones  or  heating  constraint  contacts.  The  desirable  user  imple¬ 
mentation  is  commensurate  with  the  autonomy;  namely  repeatability  and  no  a  priori 
solution  requirements.  The  next  measure  of  mission  suitability  is  the  adaptability  to 
alternate  missions.  Due  to  the  ability  to  set  up  independent  phases,  the  number  of 
waypoints  should  not  hinder  solution  performance.  Also,  since  the  knowledge  of  the 
number  of  no-fly  zone  and  heating  constraint  contacts  is  not  required  in  the  problem 
setup,  this  technique  dynamically  determines  the  optimal  solution.  Lastly,  user  im¬ 
plementation  does  not  require  extensive  training  or  excessive  intervention  to  converge 
to  a  solution.  Knowledge  of  optimal  control  theory  is  not  required,  and  the  described 
problem  setup  minimized  user  intervention  or  iterations.  Thus,  once  the  cost  and 
dynamics  are  defined,  the  phases  can  be  automatically  separated  based  on  the  given 
waypoints  and  if  a  feasible  solution  exist,  observed  convergence  has  been  extremely 
fast.  Addressing  infeasibility  is  discussed  in  Section  5.3,  while  the  contributions  of 
this  research  are  covered  in  the  next  section. 
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5. 2  Contributions 


The  contributions  arise  from  both  the  2-D  model  research  and  the  3-D  model 
research.  The  optimal  geometric  solution  is  derived  within  the  2-D  model  research. 
The  geometric  solution  is  invaluable  since  it  provides  the  optimal  solution  with  no 
integration  for  the  constant  speed  case;  and  provides  a  near-optimal  solution  for  the 
deceleration  case.  Thus,  any  autonomous  unmanned  vehicle  with  limited  comput¬ 
ing  capacity  could  benefit  greatly  from  such  a  computationally  efficient  algorithm. 
Furthermore,  any  hypersonic  vehicle  with  increased  crossrange  capability  could  ben¬ 
efit  by  using  the  geometric  solution  as  an  initial  guess  satisfying  the  waypoint  and 
no-fly  zone  constraints;  however,  the  “point-to-target”  initial  guess  in  this  research 
was  sufficient  thus  the  geometric  solution  was  not  required.  The  application  of  the 
numerical  collocation  solution  technique  to  this  highly  constrained  reentry  trajectory 
problem  is  an  outcome  of  the  3-D  model  research.  Applying  the  solution  technique 
to  the  constrained  reentry  trajectory  problem  was  one  contribution,  while  verifying 
the  optimality  was  another  equally  important  contribution.  Furthermore,  detailing  an 
efficient  user  implementation  strategy  strengthened  the  applicability  of  this  numerical 
technique.  For  example,  the  proposed  user  implementation  solution  took  minutes  on 
a  current  day  desktop  computer,  as  compared  to  hours  of  manually  coding  additional 
phase  breakpoints  and  interpreting  multiple  solution  iterations.  The  autonomy  of  the 
proposed  user  strategy  also  contributes  to  rapid  war-gaming  analysis  for  assessing 
mission  capability  and  anticipating  mission  success.  Overall,  this  research  developed 
enabling  solution  methods  and  implementation  techniques  to  enhance  Global  Strike 
mission  success. 

5.3  Recommendations  for  Future  Research 

This  research  problem  statement  intentionally  included  a  narrow  scope  of  con¬ 
straint  definitions;  motivated  by  existing  user  requests  and  compatibility  to  detailed 
analytical  techniques.  Increasing  the  model  fidelity,  broadening  the  constraint  defi¬ 
nitions,  or  increasing  the  depth  of  analysis  of  the  existing  problem,  can  open  up  a 


wide  band  of  future  research.  The  flat  Earth  model  should  be  expanded  to  a  rotating 
spherical  Earth  model  to  produce  more  operationally  representative  results.  Also,  the 
total  heating  should  be  constrained  or  verified  to  be  within  the  vehicle  design  toler¬ 
ance.  Other  topics  could  include  quantifying  the  ability  to  accommodate  additional 
waypoints  and  no-fly  zones,  or  analyzing  the  effects  of  raising  or  lowering  the  heating 
constraint.  Also,  research  could  investigate  techniques  to  counter  infeasibilities;  such 
as  creating  “soft”,  or  cost  based,  waypoints  and/or  no-fly  zones;  or  adding  the  ability 
to  increase  initial  altitude/velocity/energy  if  the  target  is  unreachable.  The  problem 
setup  in  this  research  can  easily  be  modified  to  maximize  velocity  for  a  set  terminal 
altitude,  or  maximize  terminal  energy  at  a  given  coordinate.  Altitude  could  be  added 
to  no-fly  zones  by  modeling  as  hemispheres  or  cylinders  with  a  top  altitude.  Other 
shapes  for  no-fly  zones  could  be  specified,  e.g.  ellipses  or  a  generic  shape  with  smooth 
corners,  ( x/a)p  +  ( y/b)r  =  1.  Additional  state  constraints  at  waypoint  passage  could 
be  added,  e.g.  zero  bank  angle  for  a  reconnaissance  platform.  Terminal  flight  path 
or  heading  angle  may  be  confined  to  accommodate  the  terminal  guidance  of  another 
reentry  vehicle.  Performing  a  study  of  the  number  of  nodes  as  a  function  of  com¬ 
putational  time,  versus  accuracy  and  convergence  may  also  be  of  interest.  Dynamic 
meshing  (number  of  nodes)  may  lead  to  faster  solutions;  e.g.  start  with  a  course  grid, 
then  use  that  as  the  initial  guess  to  a  finer  more  accurate  grid.  Dynamic  phase  break¬ 
points  could  be  implemented,  i.e.  adding  phases  at  identified  areas  of  interest  from 
a  previous  solution.  For  a  vehicle  with  greater  crossrange,  the  geometric  approach 
may  contribute  significantly  to  convergence  when  used  as  an  initial  guess.  Mission 
analysis  could  be  performed;  e.g.  target  footprints,  time-to-respond,  and  maximum 
accessible  intermediate  targets.  No-fly  zones  could  be  treated  more  like  threat  zones 
with  likelihood  of  damage  proportional  to  threat  proximity,  or  probability  of  survival 
could  be  calculated.  Bank  angles  in  excess  of  90  degrees  and/or  negative  angle-of- 
attack  (c(  <  0)  could  be  used  to  define  new  missions,  such  as  closer  range  with  higher 
kinetic  kill  capability.  Also,  pop-up  threats  or  real-time  waypoint  tasking  can  be  in¬ 
vestigated.  This  is  commensurate  with  the  ongoing  RTOC  research.  The  definition  of 


87 


“soft”  waypoints  could  include  achieving  the  highest  number  of  some,  yet  completely 
ignoring  others;  e.g.  maximizing  the  number  of  reconnaissance  targets  per  mission. 
The  inner-loop  guidance  could  be  implemented;  i.e.  compute  the  optimal  trajectory 
then  perform  closed-loop  tracking  of  that  trajectory.  In  conclusion,  the  above  list 
represents  yet  another  contribution  of  this  research  to  enable  a  whole  realm  of  Global 
Strike  mission  enhancements. 


Appendix  A.  Common  Aero  Vehicle  (CAV)  Aerodynamics 

A.l  CAV  Aerodynamic  Data 

The  following  is  taken  from  [69]  and  is  used  to  create  a  Mach  independent  model 
for  use  in  this  research. 


Table  A.l:  CAV-H  Aero  Data  Base 


Lift  to  Drag  Ratio  (L/D) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

2.2000 

2.5000 

3.1000 

3.5000 

3.3846 

3.2692 

3.2000 

15° 

2.5000 

2.6616 

2.9846 

3.2000 

3.0846 

2.9692 

2.9000 

20° 

2.2000 

2.3616 

2.6846 

2.9000 

2.7846 

2.6692 

2.6000 

Coefficient  of  Lift  (Cl) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

0.4500 

0.4250 

0.4000 

0.3800 

0.3700 

0.3600 

0.3500 

15° 

0.7400 

0.7000 

0.6700 

0.6300 

0.6000 

0.5700 

0.5570 

20° 

1.0500 

1.0000 

0.9500 

0.9000 

0.8500 

0.8000 

0.7800 

Coefficien 

of  Drag  (CD) 

AOA 

Mach  3.5 

Mach  5 

Mach  8 

Mach  10 

Mach  15 

Mach  20 

Mach  23 

10° 

0.2045 

0.1700 

0.1290 

0.1090 

0.1090 

0.1090 

0.1090 

15° 

0.2960 

0.2630 

0.2240 

0.1970 

0.1950 

0.1920 

0.1920 

20° 

0.4770 

0.4230 

0.3540 

0.3100 

0.3050 

0.3000 

0.3000 

CAV- hi  Aero  Reference  Area  Sref=  750  in2 
CAV- hi  Mass  m  =  2000  lbs  =  907.186  kg  (assuming  g  =  32.174  ft/s2) 


A. 2  CAV  Aerodynamic  Model 

The  coefficient  of  drag  (Co)  is  assumed  to  be  a  parabolic  function  of  the  coef¬ 
ficient  of  lift  (Cl)  [111]: 

CD=CDo  +  KC2L  (A.l) 

Taking  this  definition  of  Co  leads  to  an  expression  for  lift-over-drag,  L/D\ 

L  _CL  _  CL 
D  CD  CDo  +  KC\ 
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(A.2) 


In  order  to  find  the  Cl  that  will  produce  the  maximum  L/D,  take  the  partial  deriva¬ 
tive  of  Equation  (A. 2)  with  respect  to  Cl'- 

d(L/D)  1 _ 2KCj 

dCL  CDa  +  KCl  (CDo  +  KClf  y  '  1 

Set  Equation  (A. 3)  equal  to  zero  and  solve  for  the  optimal  Cl,  C*l\ 


CDo  +  KC2L  ~  2KC2L  =  0 


(A.4) 


Thus 


C*r  = 


C 


Do 


K 


Input  this  equation  into  Equation  (A.l)  to  find  the  corresponding  Cd,  C*d\ 


C*D  =  CDo  +  K 

C*  O/^ 

D  ~  Z(^D0 


c 


Do 


K 


Then  the  maximum  L/D  is: 


Z7* 

E  =c~ 


D 


The  normalized  coefficient  of  lift  is  defined  as: 


Q  = 


Cl 

C* 


(A.5) 


(A.6) 

(A.7) 

(A.8) 


(A.9) 


The  values  in  Table  A.l  are  plotted  in  Figure  A.l  with  the  curve  fit  using: 


Assuming  CLmi 


C*L  =  0.45 

(A.  10a) 

E*  =  3.24 

(A. 10b) 

0.9  from  Figure  A.l 

CL 

J-'max  r>  r\ 

Cd-max 

(A. 10c) 
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Figure  A.l:  CAV-H  Aerodynamic  Data  and  Model 


From  Figure  A.l  the  Mach  independence  appears  valid;  therefore,  for  this  research 
Equation  (A. 10)  is  used  to  represent  the  CAV  model.  In  Figure  A. 2  the  new  nondi- 
mensional  variable  eg  is  plotted  to  verify  the  eg  —  l  corresponds  to  the  maximum 
L/D. 


Figure  A. 2:  CAV-H  Nondimensional  Variable  Cg 
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Appendix  B.  N ondimensionalization 


Nondimensionalization  simplifies  the  equations  of  motion  and  ensures  all  states, 
paths,  and  cost  are  the  same  order  of  magnitude.  Scaling  is  critical  for  successful 
implementation  of  numerical  solution  techniques.  This  research  uses  a  flat  Earth 
model;  however,  a  radius  is  used  for  normalization  to  be  consistent  with  other  refer¬ 
ences  [41;  109;  110],  and  to  make  an  easy  transition  to  a  spherical  Earth  model  for 
future  research.  This  radius  is  the  radius  of  the  Earth  Rq  plus  the  flat  Earth  altitude. 
The  nondimensionalization  herein  starts  with  the  spherical  rotating  Earth  equations 
of  motion  and  simplifies  these  equations,  and  changes  notation,  to  arrive  at  the  flat 
Earth  equations  of  motion  used  in  this  research. 


B.  1  Spherical  Earth  3-D  Reentry  Equations  of  Motion  -  Rotating  Earth 

Assume  thrust  (T)  is  aligned  with  vehicle  velocity  vector.  The  states  are  radius 
(r)  from  the  center  of  the  Earth,  longitude  (9s),  latitude  (0s),  velocity  (Vs),  flight 
path  angle  (ys)  as  defined  in  Figure  B.l,  and  heading  (0s)  measured  from  East 
counterclockwise.  Equation  (B.l)  are  the  Kinematic  Equations  and  Force  Equations. 


es 

0S 


r  =Vs  sin  7s 

Vs  cos  7s  cos  0s 


r  cos  0s 

Vs  cos  7s  sin  0s 


(B.la) 

(B.lb) 

(B.lc) 


frp  _  J)  \ 

Vs  = - - - g  sin  7s  +  ruff  cos  0s  (cos  0s  sin  75  —  sin  0s  sin  0s  cos  7s)  (B.ld) 
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Ls  Vl 

—  cos  cr  —  g  cos  7s  H - cos  7 s 

m  r 


2 Vso;e  cos  0s  cos  0s  +  roJL  cos  0s  (cos  0s  cos  7 s  +  sin  0s  sin  0s  sin  75) 


i,s=vs 


Ls  sin  or  Vf 

- cos  7s  cos  0s  tan  0s 

m  cos  7s  r 


2 Vs^©  (sin  0s  cos  0s  tan  7 s  —  sin  0s) 


ru% 


cos  7s 


sin  0s  cos  0s  cos  0s 


(B.le) 


(B.lf) 
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Here  velocity  (Vs)  is  relative,  i.e.  velocity  is  not  inertial.  When  the  kinematic  equa¬ 
tions  are  integrated,  they  provide  the  position  of  the  vehicle  as  seen  from  the  rotating 
planet  [110].  The  Earth’s  rotation  is  a;®,  lift  is  Ls,  drag  is  D$,  bank  angle  is  a, 
and  mass  is  m.  The  subscript  S  refers  to  spherical,  thus  alleviating  confusion  with 
previously  used  variables.  The  first  term  in  Equation  (B.ld)  is  the  net  acceleration 
from  the  difference  between  the  thrust  T  and  the  drag  D$,  whereas  the  second  term 
is  the  component  of  gravity  in  the  velocity  direction.  The  terms  involving  the  Earth’s 
rotation  a;®  are  later  eliminated  with  the  zero  Earth  rotation  assumption.  The  first 
three  terms  of  Equation  (B.le)  are  the  change  in  flight  path  angle  due  to  the  vertical 
component  of  lift,  the  component  of  gravity  perpendicular  to  the  velocity  vector,  and 
the  component  of  centripetal  force  perpendicular  to  the  velocity  vector.  The  first  two 
terms  of  Equation  (B.lf)  are  the  component  of  lift  perpendicular  to  the  r-V  plane, 
and  centripetal  force  component  perpendicular  to  the  r-V  plane.  The  r-V  plane  is 
shown  in  Figure  B.2. 


Figure  B.l:  Fight  Path  Angle  [41]  Figure  B.2:  Schematic  of  the  r-V  plane  [41] 
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B.  2  Spherical  Earth  3-D  Reentry  Equations  of  Motion  -  Non- Rotating 
Earth 

The  first  step  is  to  assume  no  rotation  of  the  Earth,  ou®  =  0,  thus  simplifying 
Equations  B.l  [41;  109;  110]. 


Os 

4>s 

Vs 

is 

i>s 


Vs  sin  7 s 
Vs  cos  7 s  cos  f>s 
r  cos  (fis 

Vs  cos  75  sin  f>s 
r 


(T  -  Ds) 

m 


g  sin  75 


Ls  Vf 

—  cos  a  —  g  cos  7 5  H - cos  75 

m  r 

Lssmcr  Vf  ' 

- COS  7 5  COS  IPs  tan  05 

m  cos  7 5  r 


(B.2a) 

(B.2b) 

(B.2c) 

(B.2d) 

(B.2e) 

(B.2f) 


The  states  are  radius  from  the  center  of  the  Earth  (r),  longitude  (Os),  latitude  (05), 
velocity  (E5),  flight  path  angle  (75),  and  heading  (tfis)  measured  from  East  counter¬ 
clockwise.  The  flight  path  angle  is  defined  in  Figure  B.l.  The  radius  r  is  the  sum  of 
the  Earth  radius  R@  and  the  altitude  above  the  Earth. 


B.3  Flat  Earth  3-D  Reentry  Equations  of  Motion  -  Non-Rotating  Earth 

The  following  equations  of  motion  are  from  [114:pg.249];  however,  with  appro¬ 
priate  assumptions  they  can  also  be  derived  from  Equation  (B.2).  The  following 
assumptions  are  used  to  obtain  the  dimensional  flat  Earth  equations  of  motion: 

1.  Flight  path  angle  remains  small,  thus  cos  75  ~  1  and  sin  75  ~  75. 

2.  Thrust  is  zero,  T  =  0. 

3.  Drag  is  the  dominate  deceleration  term,  Ds  mg  sin  75.  The  term  mg  sin  75 

is  the  component  of  vehicle  weight  in  the  velocity  direction.  This  component  is 
assumed  small  compared  to  the  drag  Ds  since  the  flight  path  angle  is  small. 
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4.  Tangential  centripetal  acceleration  is  zero,  same  affect  as  being  at  the  equator, 

< h  =  0. 

5.  Change  from  initial  to  final  altitude  is  negligible  compared  to  the  radius  of  the 
Earth,  r  ^  r$ 

6.  Initial  gravity  remains  constant,  g  ~  go- 

7.  Horizontal  motion  is  only  Northing  y  and  Easting  x  from  initial  position,  thus 
±p  —  Vp  cos  dp  and  ijp  =  Vp  sin  9p. 

The  constants  r0  and  go  are  initial  altitude  and  gravity,  respectively.  The  new  vari¬ 
able  6p  within  the  assumptions  is  no  longer  longitude,  it  is  heading  angle  measured 
counterclockwise  from  East,  thus  analogous  to  ips  if  at  the  equator.  Let  hp  =  r  —  Rq> 
where  is  the  radius  of  the  Earth,  thus  hp  =  r.  Also,  the  subscript  F  refers  to  flat 
to  alleviating  confusion  with  previous  variables.  Lastly,  the  subscript  on  certain  vari¬ 
ables  are  changed  for  consistency;  Vp  =  Vs,  Lp  =  Ls ,  Dp  =  Ds,  and  7 p  =  7 5.  The 
following  are  the  dimensional  flat  Earth  equations  of  motion  with  the  assumptions 
applied: 


Xp 

=  Vp  cos  Op 

Vf 

=  Vp  sin  6p 

Hf 

=  Vf^f 

Vf 

Df 

m 

7  F 

Lp 

go 

—  cos  0 

771  Vp 

Vp 

Op 

Lf 

sin  a 

mVp 

(B.3a) 

(B.3b) 

(B.3c) 

(B.3d) 

Yl 

(B.3e) 

ro 

(B.3f) 
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B.4  Nondimensionalization 


The  dynamics,  controls,  and  path  constraints  must  be  nondimensionalized.  The 
dynamics  and  path  constraint  nondimensionalization  is  covered  in  the  following  sec¬ 
tions,  and  the  controls  a  and  q  are  already  nondimensional. 

B.4- 1  Equations  of  Motion.  The  initial  dimensional  radius  and  gravity  are  r$ 
and  go,  respectively  [41;  109].  From  [111]  h  =  0  is  defined  as  the  initial  nondimensional 
altitude;  therefore,  h  =  {rs  —  r0)/ro  or  h  =  (hF  —  hF0)/r0.  The  atmospheric  density 
p  is  assumed  to  be  a  simple  exponential,  i.e.  not  broken  up  into  multiple  layers.  The 
density  is  psi  at  the  surface  {rs  =  -R©),  and  is  po  at  the  initial  altitude  h  —  0,  thus: 

(B.4) 
(B-5) 
(B.6) 

The  density  decay  constant  is  0  —  0.14  km-1  for  this  research.  The  next  step 
is  to  nondimensionalize  Lp  and  Dp.  Let  nondimensional  V  be  defined  as  V  = 
Vp / ^g0r0  [41;  109].  This  normalization  factor  ^/goE)  represents  the  circular  or¬ 
bit  speed  at  the  initial  conditions.  From  Appendix  A  q  =  CF/Cf,  E*  =  C*L/C*D , 
C*D  =  2 CDo,  and  C*L  =  yJCDo/K  gives: 

CL  =  Q  Of 
CD  =  CDo+KC2l 

=  Cp/2  +  Kcj  (Cf)2 
=  Cb/2  +  Kcj^ 

=  C'D/2  +  cjC'n/2 
=  C'D(  1  +  cf)/2 

Q(1  +  4) 

2  E* 


(B.7a) 


(B.7b) 


p  =  psie-f>h*  =  psle~^roh+hF  o)  =  psle-^r°-R^e-0roh 
Po  =  Psie-^0  or  psi  =  poe^™ 


p  =  PoephFOe^hF  =  p0e-^hF~hFo)  =  poe~pr°h 
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Using  Equations  (B.6)  and  (B.7)  in  the  equations  for  lift  and  drag: 


LF 

DF 


-pVpCLSref 

lpoe-prohgor0V2ceClSref 

\PV2pCDSref 


:Poe 


l)T"%r„V2(PP^-Sre/ 


(B.8) 


(B.9) 


The  constant  Sref  is  vehicle  dependent  reference  area,  specified  in  Table  A.l  for  the 
CAV.  The  nondimensional  constant  B  =  (i p0roSrefCp)/(2m )  is  used  to  simplify  the 
hrst  terms  in  Equations  B.3d,  B.3e,  and  B.3f.  In  computing  B  the  variables  po, 
t0,  Sref,  and  m  must  all  be  expressed  in  consistent  units  in  order  to  cancel  out  the 
dimensionality.  The  term  Dp/m  is  an  acceleration,  and  the  term  Lp/ {mVp)  has  units 
of  1/time. 

Dp 
m 


Lf 
mVp 

(B.lOc) 


l/2p0e-^hg0r0V2C*L(l  +  c2) 
2  E*m 

BV2e~f3r°h(l  +  c|) 
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(B.lOa) 


Poe  0r°hgoroV2ceClS  ref 

2mV 

BVe-Pr0hc  [90 

V  r0 


(B.lOb) 


Substituting  the  quantities  from  Equation  (B.10)  into  Equation  (B.3)  leads  to  the 
dimensional  equations  below;  however,  they  are  now  in  terms  of  the  nondimensional 
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variables  V,  h,  and  9: 


Xp 


Vf 

hp 
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7  f 
9p 


V  y/g0r0  COS  9 

V  y/% 7)  sin  9 
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(B.lla) 

(B.llb) 

(B.llc) 

(B.lld) 

(B.lle) 

(B.llf) 


The  constants  g0  and  r0  are  used  to  nondimensionalize  Equation  (B.3),  using  length/r0, 
acceleration/ g0,  speed/^gcTo,  and  time/ y/r0/g0.  This  normalization  factor  y/gorg 
represents  the  circular  orbit  speed  at  the  initial  conditions.  The  angles  remain  in 
radians.  The  derivatives  as  a  function  of  states  and  time  must  also  be  normalized, 
e.g.  dx/dt  must  be  multiplied  by  l/r0  to  normalize  the  state  and  yJrQ/ g0  to  normalize 
time  in  the  denominator,  which  combines  to  l/A/fi,oro-  Using  a  similar  technique  for 
the  remaining  derivatives  in  Equation  (B.ll)  leads  to  the  following  nondimensional 
equations  of  motion  used  in  this  research: 


X  = 

V  cos  9 

(B.12a) 

y  = 

V  sin  9 

(B.12b) 

h  = 

V~i 

(B.12c) 

v  = 

BV2e-pr0h  ^  +  c2) 

(B.12d) 

2  E* 

7  = 

BVe-^cecosa  -  ^  +  V 

(B.12e) 

9  = 

BVe-prohcE  sin  a 

(B.12f) 

The  product  (5tq  is  dimensionless  and  is  maintained  as  two  variables  to  be  consistent 
with  previous  research  [41;  109;  110].  For  example,  for  typical  values  of  j3  =  0.14  km-1 
and  r0  =  6500  km  the  product  equals  910. 


B.f.2  Path  Constraint.  The  heating  constraint  is  nondimensionalized  using 
the  maximum  allowable  heating;  however,  it  must  be  expressed  in  terms  of  the  nondi- 
mensional  units  in  order  to  be  calculated.  Equation  36  uses  Vdim  to  generically  refer 
to  a  dimensional  variable  such  as  Vp.  Using  Equation  (B.4),  the  dimensional  heat 
rate  at  stagnation  is  [15;  41;  109;  114]: 
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(B.13) 


The  constant  kdim  is  based  on  the  heating  model,  and  rnose  is  the  vehicle  nose  radius. 
kdime-^r°-R(s)/2 

Let  K  =  f i  -  complete  the  expression  of  the  nondimensional  heating 

y^nos  eQsmax 

constraint,  as  found  in  Chapter  III  as  Equation  (37)  repeated  below: 


q,  =  Ke-^b^V3 


(B.14) 


B.f.3  2-D  Equations  of  Motion.  The  two-dimensional  equations  of  motion 
are  a  subset  of  the  dimensional  flat  Earth  equations  in  Equation  (B.ll).  Since  the 
altitude  is  constant,  h  =  0  and  7  =  0.  Also,  the  vertical  component  of  lift  must  be 
equal  to  the  vehicle  weight,  thus  Lp  =  mg/  cos  a.  The  dimensional  acceleration  Vp  is 
simply  represented  as  a^m-  The  2-D  dimensional  equations  of  motion  are  thus: 


Xp 

=  Vp  cos  9p 

(B.15a) 

Vf 

=  Vpsindp 

(B.15b) 

dp 
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Vf 
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At  constant  altitude,  g  in  Equation  (B.15c)  remains  equal  to  the  initial  gravity  defined 
as  go-  Let  the  normalized  control  u  be  defined  as: 


tan  a 

u  =  - 

tan  <7max 

Using  this  definition  for  u  and  the  same  nondimensionalization  variables  as 
tion  B.4.1,  the  nondimensional  2-D  (HCV)  equations  of  motion  are: 


x  —  V  cos  6 
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V  sin  6 
tan  <7max 
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These  are  the  equations  of  motion  defined  in  Section  3.3.5. 


(B.16) 
in  Sec- 
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Appendix  C.  Shortest  Path 

An  extension  of  [23]  is  to  include  a  secondary  point,  thus  optimizing  a  path  from 
an  initial  point,  through  an  intermediate  waypoint,  and  onto  a  final  destination.  A 
solution  to  the  two  point  problem  in  given  first,  it  will  then  be  applied  to  solve  the 
three  point  problem. 

C.l  Shortest  Path  between  Two  Points 

The  following  are  specified:  an  initial  coordinate,  [x0,2/o],  an  initial  heading, 
9q,  and  a  final  coordinate  [xf,yj].  The  time  history  within  the  initial  turn  is  from 
0  <  t  <  Ta  where  Ta  is  the  total  time  in  the  turn: 

x(t)  =  xq - sin  90  H —  sin(#0  +  ut) 

u  u 

yit)  =  2/0  H —  cos  9o - cos(90  +  ut) 

u  u 

Following  the  turn  is  a  straight  leg  lasting  Tb,  thus  time  is  Ta  <  t  <  (Ta  +  T&): 

x(t)  =  x(Ta)  +  cos  9 fit  -  Ta) 

y(t)  =  y{Ta )  +  Sin  9 fit  -  Ta) 

So, 

Xf  =  x0 - sin  9n  4 — sin  9  f  +  cos  9  fTb 

u  u 

Uf  —  2/o  H — cos  6*o - cos  9 1  +  sin  9rTb 

u  u 


(C.la) 

(C.lb) 


Multiply  Equation  (C.la)  by  sin  9f  and  subtract  Equation  (C.lb)  multiplied  by  cos  9f 


0.5 


Figure  C.l:  Two  Point  Turn  from  Initial  Coordinate  and  Heading 
Let  the  following  be  used  as  shorthand  notation: 


a  =  x  f  —  Xq  H —  sin  90 

u 

b  =  -  (yf  -  y0-  -  cos  90 
\  u 

1 

u 


c  =  — 


(C.3) 

(C.4) 

(C.5) 


This  shortens  to  a  sin  9 /  +  bcos9f  =  c,  thus  combining  the  right  hand  side  gives: 


Va2  +  b2  sin  (  9f  +  tan 


i(b 


=  c 


(C.6) 
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So  the  final  heading  (Of),  time  along  initial  arc  (Ta),  straight  leg  time  (Tb),  and  total 
time  (Tt)  are: 


Of 

Ta 


Tb 


Tt 


.-l(  c  \ 
sin  , 

Wa2  +  b2J 
Of  —  o0 

u 

a  cos  Of  —  b  sin  Of 


atan2  ( b ,  a) 


Ta  +  Tb 


(C.7) 

(C.8) 

(C.9) 

(C.10) 


C.2  Derive  Midpoint  as  Optimal  Waypoint  Position 

For  the  constant  speed  case,  the  optimal  turn  radius  orientation  is  stated  as  the 
waypoint  occurring  halfway  between  the  initiation  and  completion  of  the  turn,  see 
Figure  14(c)  in  Section  4.2.  The  geometric  derivation  of  x  =  0  is  rather  laborious; 
however,  the  following  are  two  alternative  derivations. 

Dynamic  optimization  shows  a  bang-level-bang  solution,  meaning  the  control 
is  either  zero  or  at  a  maximum  left  or  right  bank  angle.  At  constant  velocity,  a 
bang- level-bang  control  trajectory  is  seen  in  Figure  C.2.  The  initial,  waypoint,  and 
final  points  are  fixed;  [fco,  Vo],  [xw,  yw\,  and  [xf,  yf\  respectively.  The  initial  angle  0O 
and  final  angle  Of  are  free.  The  fraction  of  time  along  the  turn  that  waypoint  passage 
occurs  is  also  free,  and  must  be  determined  for  the  optimal  solution.  Let  the  waypoint 
passage  occur  at  some  fraction  (k)  between  starting  the  turn  (k  =  0)  and  exiting  the 
turn  (k  =  1);  therefore,  the  heading  angle  at  waypoint  passage  is  0o  +  k(0f  —  0o). 
For  example  in  4.2,  Figure  14(a)  is  k  =  0,  Figure  14(b)  is  k  ~  1/4,  Figure  14(c)  is 
k  =  1/2,  and  Figure  14(d)  is  k  —  1.  Let  the  time  of  flight  from  the  initial  point  to 
the  start  of  the  turn  be  Tq  and  the  time  of  flight  from  the  exit  of  the  turn  to  the  final 
point  be  Tf.  During  the  level  portions  of  flight  the  control  is  zero,  and  during  a  turn 
the  control  is  ±1  which  is  retained  generically  as  u. 
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Figure  C.2:  A  constant  speed  trajectory  from  an  initial  point  at  a  heading  90 , 

passing  though  an  intermediate  waypoint,  and  completing  the  turn  at  a  heading  6f 
to  intercept  a  final  point. 


C.2.1  Proof  Method  1  -  Equations  of  Motion.  The  following  derivation 
integrates  the  equations  of  motion  from  the  initial  point,  to  the  waypoint  depicted 
in  Figure  C.2,  then  to  the  final  point.  The  nondimensional  equations  of  motion  are 
Equation  (22)  with  V  =  0.  These  equations  can  be  integrated  from  the  initial  point 
to  the  waypoint: 
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Vw 


V2  V2 

xq  +  T0V  cos  90 - sin  6b  H - sin[d0  +  k(9t  —  do)] 
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V2  V2 

y0  +  T0V sin 90  +  — - cos 90 - - - cos[d0  +  k(9f  -  0O)] 

utan  a max  //  t  riii  crnax 
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The  second  set  of  equations  specifies  the  path  from  the  waypoint  to  final: 


xf 


Vf 


V2  V2 

xw - sin[d0  +  k(6t  —  0O)1  H - sin  Qf  +  TrV  cos  9t 

U  tan  (Trnax  U  tan  (Trnax 

(C.12a) 

V2  V2 

yw  H - cosi^o  +  k(9f  —  6*0)] - cos  6f  +  Tf  E  sin  6b 

U  tan  (Trnax  ^  tail  (Tmax 

(C.12b) 


The  leg  time  To  is  solved  from  Equation  (C.ll)  and  the  leg  time  Tf  is  solved  from 
Equation  (C.12): 
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(C.13) 


The  time  in  the  turn  is  V(9f  —  9o)/{u tan  <Tmax),  so  the  total  cost  (J)  is  the  time  from 
the  initial  point,  through  the  waypoint  in  Figure  C.2,  to  the  final  point: 


J  =  To  +  V(9f  -  90)/(utanamax )  +  Tf  (C.14) 


To  minimize  the  cost  with  respect  to  the  waypoint  passage  fraction  (fc),  the  partial 
derivative  of  J  with  respect  to  k  must  be  zero,  i.e.  dJ/dk  =  0.  Assume  9j  ^  9q  since 
that  would  imply  a  straight  line  or  no  turn.  Performing  the  partial  derivative: 
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(C.15) 
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While  in  the  turn  the  control  u,  the  speed  V,  the  maximum  bank  angle  6max ,  and  the 
term  (Of  —  0 o)  are  all  non-zero,  so  division  by  these  terms  remains  valid: 


cos  [k(0f  —  6*o)]  =  cos[(l  —  k)(0f  —  $o)] 

k(0f-0  0)  =  (1  -k)(6f-60) 

k  =  1/2 

This  shows  that  the  optimal  waypoint  passage  occurs  halfway  through  the  turn,  which 
has  the  geometric  properties  shown  in  Figure  14(c). 

C.2.2  Proof  Method  2  -  Dynamic  Optimization.  Dynamic  optimization  uses 
dynamics,  terminal,  and  intermediate  state  constraints  to  formulate  the  problem. 
With  only  one  waypoint,  and  without  considering  a  no-fly  zone,  this  is  a  simplified 
version  of  the  problem  generically  derived  in  Section  3.1.1.  The  cost  to  minimize  is 
final  time  ( tf )  and  the  heading  to  determine  for  the  optimal  solutions  is  0W,  which 
occurs  at  the  waypoint  passage  time  ( tw ).  The  nondimensional  equations  of  motion 
are  Equation  (22)  with  V  —  0.  The  terminal  and  waypoint  constraint  are  if  and  N 
respectively: 


-0  = 

1 

1 

-to 

N  = 

s, 

1 

s 

y(tw)  Uw 

The  Hamiltonian  is: 


H  =  \XV  cos  6  +  XvV  sin  6  +  A q— — Tr"mxu 

V 

The  following  are  Equations  (55)  and  (58)  without  the  velocity  dynamics  and  with 
Ps  =  0  (outside  no-fly  zones): 

Ax  =  0 

Xy  =  0 

Xe  =  XxV  sin  0  — 
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K(tf)  =  -  COS  Of/V 

\(tf)  =  -shn  0f/V 
XyV  cos  6  X e(tf )  =  0 


Equation  (64)  is  the  solution  for  the  jump  in  costates.  Since  there  is  only  one  waypoint, 
Oi  =  9W  and  tf  =  £+.  For  waypoint  passage  between  initial  and  final,  the  control  is 
constant  over  the  jump,  i.e.  u(t~)  =  u(t+).  Since  A^,  and  Xy  are  constants  (excluding 
the  jumps),  Xx(t+)  =  X x(tf)  and  A y(t+)  =  A y(tf).  At  the  start  of  the  turn  (t^),  the 
heading  is  6si  =  90.  The  jump  in  costates  Xx  and  Xy  are: 


7T  n.cr. 


^ ny 


^Oo\Uw  Usi) 

-  sin  9W  (—  cos  9f  sin  90  +  sin  9 f  cos  90] 

V  cos (9W  -  90) 

.  sin(6>/  -  dQ) 

-sm  9W— - -f- - wy 

V cos(9w  -  9 0) 

cos  9W  (  A x(t+)  sin  9si  -  A y(t+)  cos  9 ; 
cos(9w  9 


w  si ) 

cos  9W  (—  cos  9f  sin  90  +  sin  9 f  cos  90 


V  cos (9W  -  90) 
sin(6>/  -  6*q  ) 


sm( Uf  -  V0) 
cos  9W  .  . 

V  cos (9W  -  9 o) 


(C.17a) 


(C.17b) 


Therefore  the  costates  behind  the  jump  at  tw,  going  backwards  in  time,  are: 


Xx{tw ) 
Xx(tw) 


cos  9f  .  sin(6b  —  90) 
~V  SmS“Vcos(0„,-0o) 


Substitute  these  values  into  the  Hamiltonian  at  the  start  of  the  turn  (t~ ) ,  where 
=  0,  xitjj)  =  V  cos  0O,  and  =  V  then  set  the  Hamiltonian  equal  to 
— 1.  The  Hamiltonian  and  the  states  are  continuous  so  the  plus  and  minus  superscript 
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notation  is  unnecessary: 


tju  \  f  a  ■  a  sin(fy  ~°o)\  Q 

H(tai)  =  (-  cose,  -  sin8.cos(^_w  J  COS«0 

+  (-sinfy  +  cos0w  sin^o  —  — 1 

V  cos(6w-d0)J 

_  Q  \ 

—  (cos  Of  cos  0O  +  sin  9f  sin  0O)  —  (sin  9W  cos  0O  —  cos  9W  sin  90 ) - — - — 

cos(6llu  -  0O) 

cos (9f  —  0O)  cos(0„,  —  0O)  +  sin(0j  —  0O)  sin(0^  —  0O)  =  cos(0^  —  0O) 


-1 


cos (0/  -  0O  -  9W  +  0O)  =  cos^  -  0O) 


9 


W 


9f  -  9W  =  9W  -  0O 

e^  =  «o  +  ^(e/-eo) 


(C.18) 


Since  9W  is  of  the  form  90+k(9f— 0O),  the  optimal  solution  is  again  shown  to  be  k  =  1/2. 
This  solution  for  9W  can  be  plugged  back  into  the  jump  equation,  Equation  (C.17), 
to  solve  for  7inx  and  7Tny.  With  a  solution  for  the  jumps,  the  value  of  the  costates  at 
t~  can  be  computed.  For  one  waypoint  this  provides  the  value  of  the  initial  costates 
at  t0 ;  A x(t0)  =  -  cos 90/V,  Xy{t0)  =  -  sin 90/V,  and  A e(t0)  =  0. 
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Appendix  D.  Midpoint  Algorithm 


1  $ 

Up  y  '■ 


Figure  D.l:  Nomenclature  for  Iterative  Solution 


The  data  to  be  given  are  initial  position  [x0,  yo],  intermediate  waypoint  [xw,  yw\, 
final  position  [xf,  yf ],  radius  (r),  and  direction  of  turn  ( u ).  Control  ( u )  positive  im¬ 
plies  increasing  heading  angle  (6)  thus  counterclockwise  rotation.  Heading  is  measured 
counterclockwise  from  the  x-axis.  The  given  points  will  be  normalized  by  the  radius 
so  radius  (r)  is  assumed  1.  Also,  the  given  geometry  is  rotated  such  that  the  x-axis 
is  aligned  along  the  initial  and  final  points.  Lastly,  the  problem  is  translated  such 
that  the  initial  point  is  [0,  0].  This  can  all  be  reversed  following  the  computations  to 
recuperate  the  original  configuration. 

In  addition  to  the  given  values,  other  computed  values  are  fixed.  The  angles 
within  the  triangle  between  initial  and  the  waypoint  ( 6$ )  and  final  and  the  waypoint 
( Of ),  and  the  distances  from  initial  to  waypoint  (do),  waypoint  to  final  ( df ),  and 
initial  to  final  (do/)-  The  size  of  A 90  will  be  adjusted  until  the  optimal  orientation  is 
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achieved.  The  iteration  loop  is  initiated  with  a  guess.  A  good  initial  guess  is  described 
at  the  end  of  the  section. 

A  90  =  guess  >  0 

First,  the  relationship  between  the  change  in  A9q  and  the  corresponding  value  of  A9wo 
must  be  derived.  Figure  D.2  shows  the  initial  leg;  however,  it  is  translated  to  start 
at  [0,  0]  and  rotated  to  be  aligned  with  the  x-axis.  The  location  of  the  center  of  the 
circle  rotated  A9w0  about  the  waypoint  is: 


xc 

do  -  cos(|  -  A9w0) 

d0  -  cos(A 9w0  -  f ) 

d0  —  sin  A9w0 

iJc 

—  sin(|  -  A 9wo) 

sin(A6U)  -  f ) 

—  cos  A^o 

(D.l) 


Originating  from  the  center  of  the  circle  [xc,  yc],  at  a  perpendicular  distance  of  radius 
r  =  1,  there  exist  an  angle  9  such  that  the  line  passes  through  the  initial  point  [0,  0]. 
The  normal  form  of  the  straight  line  equation  for  any  position  [x,  y\  is: 

(x  —  xc )  cos  9  +  (y  —  yc )  sin  9  =  1  (D.2) 

In  Figure  D.2  the  right  triangle  containing  A 90  has  the  complementary  angle  equal  to 
|  —  A9q.  Since,  it  is  on  the  x-axis  which  is  parallel  to  the  horizontal  line  at  the  center 
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of  the  circle  then  it  is  also  equal  to  n  —  9.  Combining  these  relationships  provides: 


|-A  d0  =  7 r-d  (D.3) 

0  =  Ad0  +  ^  (D.4) 

Now  substituting  [x,  y]  =  [0,  0],  Equations  (D.l)  and  (D.4)  into  Equation  (D.2) 
produces: 


(0  -  d0  +  sin  A9w0)  cos(A d0  +  ^-)  +  (0  +  cos  A 9w0)  sin(Ad0  +  ^)  =  1  (D.5) 

(d0  —  sin  A9w0)  sin  A 90  +  cos  A9w0  cos  A 90  =  1  (D.6) 

Equation  (D.6)  establishes  the  relationship  between  A 90  and  A9wq.  Since  the  itera¬ 
tions  will  start  with  A d0  the  expression  for  A 9w0  in  terms  of  A 90  must  be  found: 


d0  sin  A 90  +  cos  A 9w0  cos  Ad0  —  sin  A 9w0  sin  A 90 

cos(Adu,o  A  Ad0) 
Ad„,o  =  cos_1(l  —  do  sin  Ad0)  —  Ad0 


1 

1  —  d0  sin  Ad0 

(D.7) 


In  Figure  D.l,  at  the  waypoint,  the  relationship  between  A9wq  and  A 9wf  is  defined 
as: 


A9wq  +  A9wf  +  6  —  7r 

S  =  TT  —  9q  —  9  f 
A9w0  +  A9wf  +  7T  -  d0  -  9f  =  n 
A  9wf  =  90  +  9f  —  A9w0 


(D.8) 


Now  the  process  needs  to  be  reversed  for  the  final  leg.  The  same  setup  as  in  Figure 
D.2  can  be  used,  except  all  of  the  0  subscripts  become  /  subscripts.  The  terms  will 
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be  regrouped  differently  in  order  to  get  A  9j  in  terms  of  the  computed  A  9wf 


XCf 

df  —  sin  A 9wf 

—  cos  A 9wf 

xCfsmA9f  —  yCfcosA9f  =  1 


dc  =  yJx2Cf  +y , 


dc  sin  I  A 9f  +  tan 


-l 


Ucf 


x, 


cf 


sin  yA9f  —  tan  1  ^  — 


=  1 
1 

dr 


(D.9) 


Using  the  Law  of  Sines,  the  exterior  triangle  can  be  computed  from  the  computed 
interior  angles  and  known  side  length  do/: 


A  —  9  o  +  A9q 

B  =  9f  +  A9f 

C  =  t r-A-B 
b  _  d0f 


sin  B  sin  C 

i  dof. 

b  =  — —  sm  B 
sin  C 


xp 

bcosA 

Up 

b  sin  A 

The  optimal  solution  occurs  when  the  line  from  the  center  of  the  circle  [xc,  yc\  to  the 
peak  of  the  triangle  [xp,  yp]  passes  through  the  waypoint  [xwi  yw\.  The  center  of  the 
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circle  in  Equation  (D.l)  must  be  rotated  to  realign  with  the  initial  leg: 


xc 

!)c 


0r  =  tan 


cos  60  —  sin  6  0 
sin  60  cos  60 

—i  (  Up  He 


~i 


6b. .  =  tan 


error  =  6b„  -  6b 


Uw  Vc 


Xn, 


Xr 


The  iteration  continues  with  a  new  value  for  A6b 


xc 

Vc 


A0o  =  A  6*o  +  error 


(D.10) 


(Dll) 


(0.12) 


This  iterative  loop  was  compared  to  satisfying  Equation  (D.ll)  using  Matlab®’s 
f  solve.  Accuracy  was  driven  to  le-10  radians.  Although  the  loop  took  20  to  30 
iterations  compared  to  only  2  to  4  with  f  solve,  the  iteration  loop  took  orders  of 
magnitude  less  epu  time.  The  other  advantage  is  f  solve  sometime  took  excessively 
large  jumps  in  A6q  which  created  infeasible  constructions,  i.e.  trigonometric  function 
errors.  The  final  step  is  to  convert  the  iteration  results  to  vehicle  headings.  For  this 
configuration: 


A 

-B 

Of  +  6*o 
2 

Equation  (D.9)  can  be  used  to  compute  an  initial  guess  for  A60.  A  good  guess  is  that 
the  radius  will  be  aligned  with  the  bisect  of  angle  5.  In  many  cases  this  initial  guess  is 
within  a  degree  of  the  final  answer;  therefore,  a  non-iterative,  near  optimal,  solution 


0  o  = 

df  = 
0 = 
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exists  to  minimize  computational  time: 


a  e. 


wOR 


Xr. 


y  c 


dr 

t-guess 

^^Oguess 


7 r  5  Oq  +  Of 
2~  2~  2 


dn  —  sin  Ad , 


wOs 


cos  A0wOs 


xt 


sin 


-1 


+  yl 


dr 


+  tan 


-i 


xr 
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Appendix  E.  Pseudo  spectral  Method 


The  purpose  of  pseudospectral  methods  is  to  approximate  the  continuous  solu¬ 
tion  to  a  set  of  differential  equations  using  polynomial  interpolation  through  discrete 
points  or  nodes.  The  motivation  is  to  avoid  sequential  integration  which  can  lead  to 
divergence  and  may  prohibit  determining  a  solution.  Legendre  or  Chebyshev  polyno¬ 
mials,  as  shown  in  Figure  E.l  [31],  may  be  used  to  satisfy  the  differential  equations 
at  a  discrete  number  of  nodes  N .  Computing  the  solution  at  the  nodes  is  also  termed 
collocation  [11;  31;  50;  99].  These  collocation  techniques  satisfy  all  the  nodes  simulta¬ 
neously,  thus  avoiding  the  pitfalls  of  integration,  especially  the  forward  and  backward 
integration  within  the  shooting  method.  The  discretization  and  node  placement  is 
fundamental  to  the  various  pseudospectral  techniques  [4;  33;  80;  115;  116].  Equally 
space  nodes  is  the  simplest  arrangement;  however,  for  polynomial  interpolation  equi- 
spacing  will  lead  to  large  errors  as  the  number  of  nodes  is  increased.  This  large  error, 
especially  at  the  endpoints,  is  called  the  Runge  Phenomenon  as  seen  in  Figure  E.2(a). 


Figure  E.l:  Comparison  of  Legendre  and  Chebyshev  Polynomials  [31] 
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error  =5.7818 

max 


(a)  Evenly  Spaced  Nodes 

error  =0.0175 

max 


(b)  Chebyshev  Nodes 

Figure  E.2:  Runge  Phenomenon  is  Avoided  Using  Chebyshev  Points 


The  avoid  these  errors,  the  density  of  the  nodes  at  the  endpoints  is  increased 
using  Chebyshev  points.  The  Chebyshev  point  placement  in  Equation  (E.l)  is  shown 
in  Figure  E.3. 

Xj  =  cos(jn/N),  j  =  0, 1, . . . ,  N  (E.l) 

Using  Chebyshev  points,  the  polynomial  interpolation  accuracy  improves  with  an 
increase  in  the  number  of  nodes  as  seen  in  Figure  E.2(b).  An  important  observation 
is  that  the  node  placement  is  more  sparse  at  the  middle;  therefore,  simply  increasing 
the  number  of  nodes  may  not  produce  more  accurate  results  if  the  region  of  interest  is 
near  the  midpoint.  Software  packages  such  as  DIDO  and  GPOCS  allow  for  regions  of 
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more  dense  node  placement.  This  is  implemented  using  “knots”  in  DIDO  or  “phases” 
in  GPOCS.  By  specifying  these  intermediate  events,  or  boundary  conditions,  a  new 
endpoint  is  created;  thus,  adding  nodes  and  reducing  errors. 


Figure  E.3:  Chebyshev  Points  Projected  onto  the  x-axis  [99] 

Increasing  the  number  of  nodes  does  improve  accuracy;  however,  it  may  also 
extend  the  time  to  compute  the  solution.  Also,  it  may  be  more  efficient  to  compute 
a  solution  with  a  reduced  number  of  nodes  first,  then  use  that  solution  as  an  initial 
guess  for  a  subsequent  software  run  with  additional  nodes.  This  dynamic  meshing 
technique  is  mentioned  as  potential  additional  research  in  Section  V. 

The  pseudospectral  method  requires  that  the  differential  equation  be  exactly 
satisfied  at  the  “collocation”  points  [ll:pg.l2].  This  may  be  accomplished  using 
Chebyshev  differential  matrices  [99:pg.54],  An  example  from  [99:pg.54]  is  recreated 
and  presented  in  Figure  E.4.  These  methods  use  state  values  at  discrete  nodes  to 
construct  differential  matrices  [31],  to  then  compute  the  derivative  as  shown  in  Fig¬ 
ure  E.4.  The  derivative  matrix  used  to  compute  the  results  in  Figure  E.4  is  computed 
using  N  =  7  and  IV  =  20.  This  demonstrates  the  initial  accuracy  possible  even  with 
a  lower  number  of  nodes,  and  shows  how  the  accuracy  increasing  with  an  increased 
number  of  nodes. 

The  nodes  included  in  the  solution  to  the  differential  equations  may  differ  be¬ 
tween  software  packages.  DIDO  uses  the  Legendre  pseudospectral  method  which 
is  based  on  interpolating  functions  on  Legendre-Gauss-Lobatto  (LGL)  quadrature 
nodes  [78].  This  technique  does  include  the  boundary  conditions  whereas  the  GPOCS 
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Gauss  pseudospectral  technique  using  Legendre-Gauss  nodes  [4]  specifically  does  not 
use  the  endpoints.  This  variation  in  solution  technique  also  varies  the  discrete  KKT 
multiplier  mapping  to  the  continuous  costates.  The  Covector  Mapping  Principle  [27] 
is  used  within  DIDO;  however,  the  Costate  Mapping  Theorem  [4]  is  used  within 
GPOCS.  The  placement  of  this  mapping  is  shown  in  Figure  E.5. 


Figure  E.4:  Derivative  Computed  Using  Chebyshev  Differentiation  Matrices 


Problem  B  x 


convergence 


approximation 
(indirect  method) 


Problem  B 

i ,  Covector 
Mapping 
Theorem 


convergence 

Problem  B  r  ~  Problem  B  N 

approximation 
(direct  method) 


(a)  DIDO  [27] 


(b)  GPOCS  [4] 


Figure  E.5:  Computation  of  the  Continuous  Lagrange  Multipliers 
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Appendix  F.  Software  Implementation  Strategies 

The  software  coding  techniques  used  to  compute  the  numerical  results  are  pre¬ 
sented  in  this  appendix.  This  is  included  to  help  any  follow-on  research  duplicate  the 
results  presented  herein;  however,  it  not  intended  to  replace  or  replicate  the  software’s 
user’s  manual  [72], 

F.l  Capture  Jumps  in  Costates 

The  waypoints  and  final  altitude  are  specified  as  matching  lower  and  upper 
bounds  on  the  initial  and  terminal  state  bounds  for  the  applicable  phase,  i.e.  they  are 
not  computed  as  separate  events.  The  no-fly  zones  are  specified  as  S  in  Equation  (26); 
however,  to  force  the  software  to  find  the  exact  location  of  the  discontinuities  the 
definitions  for  S'  =  0  and  S''1)  =  0  are  specified  to  be  satisfied  at  the  end  of  the 
applicable  phases,  i.e.  these  are  event  conditions  as  a  function  of  the  terminal  states 
for  that  phase.  If  higher  accuracy  is  required  to  force  a  sharp  discontinuity,  instead  of 
one  that  has  a  partially  continuous  corner,  the  S ^  function  is  multiplied  by  a  scaling 
factor  such  as  100.  The  heating  constraint  Q  —  0  and  Qx(q  =  0.8)  <  Q1  <  Q1(q  = 
1.2)  are  set  as  events  to  force  a  discontinuity  at  the  heating  constraint  contact.  In 
looking  at  Equation  (43)  it  can  be  seen  that  it  is  a  function  of  the  control,  which  is 
unknown;  hence,  the  band  selection  of  0.8  <  q  <  1.2.  Again,  these  were  tricks  to 
force  and  exact  solution,  final  implementation  is  much  simpler  by  merely  setting  the 
event  7  =  0.  This  doesn’t  force  a  jump  as  accurate  as  above,  but  it  does  provide 
sufficient  nodes  to  properly  characterize  the  jump  to  converge  on  the  correct  solution. 

F.2  Control  Multipliers 

The  simplest  method  for  control  limitations  is  to  specify  the  upper  and  lower 
bounds,  per  phase.  If,  however,  the  multipliers  associated  with  an  active  control 
constraint  are  of  interest,  the  control  can  be  input  as  a  path  constraint.  The  control 
multiplier  is  then  contained  in  the  output  path  multipliers.  This  is  the  method  used 
to  compare  the  numerical  and  analytical  results  within  this  research. 
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F.3  Events  and  Cost 


There  were  a  few  techniques  that  provided  better  results;  however,  I  did  not  run 
any  formal  tests.  First,  the  minimum  time  cost  performed  better  when  it  was  specified 
as  the  terminal  time  of  the  final  phase.  This  seemed  to  work  better  than  having  an 
integrand  cost  of  1.  Also,  be  sure  to  use  the  intended  variable  since  initial,  terminal, 
and  intermediate  values  are  sometimes  available  in  the  same  function.  Typically, 
for  events  and  linkages  the  intended  variable  is  initial  or  terminal.  Inserting  control 
limitations  and  state  constraints  worked  best  when  inserted  directly  as  part  of  the 
state  or  control  bounds;  therefore,  use  the  bounds  instead  of  a  separate  path  or  event 
if  possible. 

F.4  Singularities  and  Results 

Always  ensure  any  singularities  are  avoided;  therefore,  set  state  or  control 
bounds  to  avoid  undefined  trigonometric  functions  or  a  divide  by  zero.  Integrat¬ 
ing  the  dynamics  is  a  beneficial  way  to  ensure  the  answer  is  feasible.  Also,  ensure  the 
maximum  step  size  of  the  integrator  is  low  enough  to  avoid  divergence.  The  mag¬ 
nitude  of  the  costates  and  the  value  of  the  Hamiltonian  are  valuable  in  identifying 
an  optimal  solution.  If  the  costate  have  huge  magnitudes,  or  the  Hamiltonian  does 
not  remain  constant,  a  non-optimal  solution  is  likely.  The  Hamiltonian  will  only  be 
constant  if  it  is  not  explicitly  a  function  of  time. 

F.5  Nodes  and  Iterations 

It  may  be  best  to  start  a  problem  with  fewer  nodes  and  a  low  maximum  number 
of  major  iteration  steps.  This  technique  allows  for  an  initial  look  at  the  results  to 
verify  it  is  approaching  an  expected  or  reasonable  solution.  It  also  can  be  used  to 
identify  contact  with  a  path  constraint.  Such  contact  may  require  another  phase 
breakpoint  to  increase  the  number  of  nodes  at  that  area  of  interest.  This  may  also  be 
a  detriment,  too  few  nodes  with  high  optimality  tolerance  and/or  feasibility  tolerance 
can  produce  an  “infinite”  loop.  The  controls  may  be  so  spaced  out  that  a  small 
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change  is  held  for  such  a  long  time  that  convergence  is  not  possible,  or  at  least  is  very 
time  consuming.  Typically  this  behavior  can  be  identified  by  monitoring  the  merit 
function  and  looking  for  bounded  oscillations. 

F.6  Breakwell  Example 

The  Breakwell  problem  is  presented  in  two  publications,  but  with  different 
costate  histories  [14;  78].  The  results  in  [78]  are  reproduced  using  GPOCS  and  are 
plotted  in  Figure  F.l.  This  difference  with  respect  to  the  research  herein  is  addressed 
in  Sections  4.3  and  4.6.  This  example  demonstrates  how  to  produce  the  same  results 
using  the  analytical  method  and  numerical  method.  The  basis  of  confusion  is  the 
analytical  method  [14;  44]  adds  interior-point  constraints  to  provide  more  criteria 
for  a  solution.  This  adds  intuition  into  the  solution;  however,  it  is  not  required  for 
the  numerical  results.  Although  these  additional  steps  are  not  required,  they  will  be 
implemented  in  the  numerical  method  to  force  comparable  results.  This  comparison 
will  verify  that  even  though  these  methods  produce  different  costates,  they  are  both 
correct,  and  both  lead  to  the  same  optimal  control  and  state  history. 

F.6. 1  Analytical  Derivation.  The  differences  in  costates  occur  when  adjoin¬ 
ing  the  path  constraint  S,  or  one  of  its  time  derivatives,  to  the  Hamiltonian.  For  the 
Breakwell  problem  the  constraint  is  x  <  £,  thus  the  path  constraint  is: 

S(x(t),t)  =  [x(t)  —  Z\  <  0  (F.l) 

For  the  analytical  method,  the  next  step  is  to  define  the  interior-point  constraint  that 
is  applied  upon  entering  the  constraint  boundary  at  time  t\.  This  is  derived  by  taking 
the  time  derivative  of  S  until  the  control  a  appears.  For  this  problem  a  —  v  —  x; 
therefore,  the  second  time  derivative  contains  the  control,  i.e.  S(-2]  =  a.  Thus,  the 
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lower  time  derivatives  become  the  interior-point  constraint: 


N(x.(ti),ti) 


5 

1 

I 

sw 

v(h) 

(F.2) 


From  the  definition  of  N  there  will  be  a  jump  in  Xx  and  A„,  as  shown  in  Figure  F.l(b) 
for  the  analytical  case.  Typically,  GPOCS  does  not  compute  the  costates  in  the 
manner  just  described  since  this  is  more  derivation  than  simply  programming  the 
path  constraint  S. 


F.6.2  GPOCS  Implementation.  The  GPOCS  software  package  allows  the 
problem  to  be  broken  up  into  phases.  To  force  GPOCS  to  replicate  the  Bryson 
derivation,  the  problem  is  split  into  three  phases.  Since  the  Bryson  problem  setup  is 
going  to  be  three  phases,  the  initial  comparison  is  also  broken  into  three  phases,  as 
shown  in  Figure  F.l.  This  is  not  required  for  a  solution,  but  it  does  provide  a  better 
comparison  of  the  computed  cost. 

Figure  F.l  uses  S  as  the  adjoined  path  constraint.  For  this  particular  numerical 
solution,  the  problem  setup  also  requires  an  interior-point  constraint  of  S  =  0  at  the 
end  of  phases  1  and  2.  Thus,  since  S  is  only  a  function  of  x,  there  is  only  a  jump  in  the 
A*  costate.  The  next  step  is  to  replicate  the  Bryson  results  derived  in  Section  F.6.1. 

Figure  F.2  uses  S''2'  as  the  adjoined  path  constraint.  From  the  previous  deriva¬ 
tion  there  is  also  an  interior-point  constraint  N,  from  Equation  (F.2),  at  the  end  of 
phase  1,  i.e.  entering  the  path  constraint  boundary.  The  linkage  between  phases  is 
programmed  to  force  state  continuity,  but  the  costates  and  control  may  be  discontinu¬ 
ous.  The  path  constraint  for  the  second  phase  is  S ^  =  a  =  0.  Since  the  interior-point 
constraint  N  is  a  function  of  x  and  v,  there  is  a  jump  in  A^  and  Xv,  as  shown  in  Fig¬ 
ure  F.2(b).  This  GPOCS  problem  setup  replicates  the  analytical  Bryson  derivation; 
therefore,  the  GPOCS  and  Bryson  solutions  now  match,  as  seen  in  Figure  F.2(b). 
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Cost  D  =  4.4444444,  Cost  =  4.4444445 

Bryson  GPOCS 


Figure  F.l:  Breakwell  Problem:  Differences  in  Bryson  and  GPOCS  Results 
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Cost  D  =  4.4444444,  Cost  =  4.4444445 

Bryson  GPOCS 


Figure  F.2: 


Breakwell  Problem:  Matching  Bryson  and  GPOCS  Results 


F.6.3  Matlab®  Code.  The  following  code  contains  the  GPOCS  solution, 
the  Bryson  and  Ho  solution,  and  plotting  commands.  The  GPOCS  main  program  is 
Breakwell_Problem.  The  main  program  calls  the  cost  function  (Cost),  the  equations 
of  motion  (Dynamics),  the  path  constraint  for  the  second  of  three  phases  (Path),  and 
the  linkage  between  multiple  phases  (Link).  The  support  functions  are  provided  as 
function  handles  using  the  @  call,  thus  allowing  all  of  these  functions  to  be  contained  in 
a  single  hie.  For  comparison,  the  analytical  Bryson  solution  is  computed  in  Bryson_Ho. 
The  plotting  routine,  Plot_Results,  is  provided  to  demonstrate  how  to  extract  the 
cost,  time,  states,  costates,  Hamiltonian,  and  path  constraint  multipliers  from  the 
GPOCS  results  structure. 


function  Br eakwell_Pr oblem  7.  Multiple  Methods,  Multiple  Phases 
7. l  written  by:  Maj  Tim  Jorris,  AFIT/ENY,  DSY-07S,  Grad:  Sep  2007 
clear  global  ,  global  CONST  7.  Remove  previous  global  values 
7,7,  Constants 

CONST. ell  =  1/10;  7,  maximium  ,  x(t)  <=  ell 


7.7.  Method 

CONST . bryson.method  =  true; 

7.7.  Functions 
problem. FUNCS . cost 
problem . FUNCS  .  ode 
problem . FUNCS . path 
problem . FUNCS  .  link 

7.7.  Number  of  Phases  Required 

N  =  36;  7,  Total  Number  of  Nodes 

if  CONST  .  ell  >  1/4  7.  Unconstrained  Problem 


=  OCost  ; 

=  QDynamics 
=  OPath  ; 

=  OLink  ; 


7.  true  or  false 

7.  Cost 

7.  Equations  of  Motion 
7.  Path  Constraint 
7.  Linkage  between  Phases 


7.7.  Nodes 

N=N ;  Np=length (N) ;  CONST. Np=Np; 
for  i=l:Np,  ph ( i ). nodes =N ( i ) ;  end 

7. 7.  Time 

tf  =  l;  tbreak=[0  ;  tf]  ;  dt=[0  0]  ;  7.  breakpoints  ,  +/-  amount 

problem  .  FUNCS  =  rmf  ield  (problem  .  FUNCS  path ’)  ;  7.  not  required 
elseif  CONST,  ell  >=  1/6  7.  IP  Constraint  between  Phases  1  and  2 


7.7.  Nodes 
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Ni  =  round (N/2)  ;  N=  [Ni , Ni]  ;  Np  =  length (N)  ;  CONST. Np  =  Np; 
for  i=l:Np,  ph ( i ). nodes =N ( i ) ;  end 
ll  Time 

tf=l;  tbreak=[0  0.5  t f ]  ;  dt=[0  0  0 ]  ; %  breakpoints,  +/-  amount 

problem . FUNCS  =  rmf ield (problem . FUNCS  path.  ’ )  ;  %  not  required 
else  1  3  Phases ,  Path  Constraint  for  Phase  2 

•/,'/.  Nodes 

Ni  =  round (N/3)  ;  N=  [Ni , Ni , Ni ]  ;  Np= length (N)  ;  CONST. Np  =  Np; 
for  i=l:Np,  ph ( i ). nodes =N ( i ) ;  end 
'/.  it  Time 

tf  =1 ;  tbreak=  [0  3*C0NST.ell  1 -3*  CONST . ell  tf ]  ;  dt=[0  0  0  0]; 

end 

7, 1  Time  of  Phase  Breakpoints 
for  i = 1 : Np 

ph(i) .time.min=[tbreak(i)-dt(i)  tbreak(i+l)-dt(i+l)] ; 
ph(i) .time. max =[tbreak(i)+dt(i)  tbreak(i+l)+dt(i+l)] ; 

end 

ii  States 

x0=0;  v0=  1;  xf=0;  vf=-l;  xnames ={ ’ x ’ , ’ v ’ } ;  unames={,a'}; 

for  i=l:Np,  ph(i) . states . names=xnames ; 
ph(i)  . states .min  =  -2*ones (2 ,3)  ; 
ph ( i ) . st at e s . max =  2*ones(2,3); 
if  i==l 

’/,  fixed  initial 

ph(i) . states. min ( : ,l)=[x0;v0] ; 
ph(i) . states. m ax (: ,l)=[x0;v0] ; 

‘/,  fixed  phase  endpoint  (if  multiple  phases) 
if  Np  >=  2 

ph(i)  .  states. min(:  ,end)  =  [C0NST.ell;0]  ; 
ph(i)  .  states.max(:  ,end)  =  [C0NST.ell;0]  ; 

end 

end 

if  i  =  =  2  &  i~  =  Np  it  only  if  multiple  phases 
it  fixed  phase  initial 

ph(i) . states .min ( : ,l)=[C0NST.ell;0] ; 
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ph  (i) . states . max ( : ,l)=[C0NST.ell;0] ; 
if  ~ CONST . bryson_method 
%  fixed  phase  endpoint 

ph(i) . states. min(: , end)=[C0NST.ell;0] ; 
ph(i) . states.max(: , end)=[C0NST.ell;0] ; 

end 

end 

if  i==Np 

if  ~ CONST . bryson_method  &  Np  >=  2 

/  fixed  phase  initial  (if  multiple  phases) 
ph(i) . states. min(: ,l)=[C0NST.ell;0] ; 
ph(i) .states.max(: ,l)=[C0NST.ell;0] ; 

end 

l  fixed  terminal 

ph(i) . states. min(: ,end)=[xf ;vf] ; 
ph(i) .states. max (: ,end)=[xf;vf] ; 

end 

end 

“/, “/.  Controls 

for  i=l:Np,  ph ( i ). controls . names =unames ; 

ph (i) . controls . min=-10 ;  ph (i) . controls . max=  10; 

end 

11  Path 

if  i sf ield ( problem . FUNCS , ’ path ’ ) 

ph(2)  .path. min  =  [0]  ;  ph(2)  . path . max  =[0]  ; 

end 

11  Guess 

if  sum(N)  <=  36  1  Initial  low  number  of  nodes  for  estimated  sol’n 
t_guess  =  tbreak (  [1 ; end] )  ;  x_g=[  xO  vO ;  xf  vf  ]; 

for  i=l : Np 

guess ( i ). t ime  =[tbreak(i)  ; tbreak ( i +1 ) ] ; 

guess ( i ). control s = int erp 1 ( t.guess ,  x_g  , guess ( i ). t ime ) ; 

guess ( i )  . st at e s  =interpl (t.guess  ,  x_g  , guess ( i ). t ime )  ; 

end 

else  1  Increased  number  of  nodes  for  higher  accuracy 
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Prev  =  load  (  [mf  ilename  ,  ’  .Results  ’  ]  )  ;  1  Previous  estimated  sol’n 

for  i=l : Np 

guess(i) .time  =  Prev . results . solution (i) . time ; 

guess ( i ). control s  =  Pr ev . result s . solut i on ( i ). controls ; 
guess  (  i  )  .  st  at  e  s  =  Pr  ev  .  result  s  .  solut  i  on  (  i  ).  st  at  es  ; 

end 

end 

11  Linkage 

for  i=l:(Np-l)  1  Always  one  less  link  than  phases 
1  Want  to  enforce  state  continuity 
1 ink age  s ( i )  . left . phas  e  _ number  =  i ; 
linkages (i) . right . p h as e_ number =i+l ; 

1 ink age  s (i).min=0*ph(l).states.min(:,l);  1  gets  right  size 
1 ink age  s (i)  ,max=linkages(i)  .min; 
problem . linkages=linkages ; 

end 

11  Problem 

problem . phases  =  ph  ;  problem . guess  =guess;  problem . name  =  mf ilename  ; 
problem.  independent_variable= ’ increasing'  ; 

problem . Derivatives=' automatic ’ ;  problem. autoscale=’off ’ ; 
problem . FeasTol =1 e -6 ;  pr oblem . OptTol = 1 e -6 ;  PrintOutput=true ; 

11  Run 

results=gpocs ( ’ snopt7 ’ , problem , PrintOutput) ; 

save ( [mf ilename , ’ .Results ' ] )  1  Used  as  estimate  for  later  runs 

Plot_Results (results) 

11  Cost  Function 

function  [phi,  L]  =  Cost(sol) 

phi  =  0;  L  =  1/2* sol . controls . “2 ; 

11  Equations  of  Motion 
function  Xdot  =  Dynamics ( sol ) 

dx  =  sol . states ( :  ,2)  ;  dv  =  sol . controls (:  ,1)  ;  Xdot  =  [dx,dv]; 

11  Path  Constraint 
function  pth=Path ( sol ) 

global  CONST  1  phase.number  is  not  necessary 
if  sol . phase.number ==2  &  CONST . bryson.method 
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pth=sol . controls ; 

else  •/.  This  function  is  only  called  for  Phase  2 

pth=sol . states ( : , 1) -CONST . ell ; 

end 

ll  Linkage 

function  link=Link (breakpoint ) 

xminus=breakpoint . right . state ;  xplus=breakpoint . left . state ; 
link  =  xminus  (1  :  end) -xplus  (1  :  end)  ;  °/0  continuous  states 
7.7.  End  of  GPOCS  Code 

ll  Bryson  and  Ho  Solution,  pg .  121-123 

f unct ion  [results] =Bryson_Ho (ell) 

N=201;  %  Define  number  of  time  steps 
if  ell  >  1/4,  t=linspace (0 , 1 , N) ; 
elseif  ell  >=  1/6 

n  =  round(N/2)  ;  7.  ensures  1/2  is  included 
t=[linspace(0,l/2,n)  ’  ;  linspace(l/2+eps,l,n)  ’  ]  ’  ; 

else 

n  =  round (N/3)  ;  t=[linspace(0,3*ell,n)  1  ;  .  .  . 

linspace(3*ell  +  eps,l-3*ell,n)  ’;linspace(l-3*ell+eps,l,n)  1  ]  1 

end 

7.  Define  states  and  costates  ,  p.121  Bryson  and  Ho 
if  ell  >  1/4 

u  =  0*t  ;  x  =  u;  v=u;  lamx  =  u;  lamv  =  u;  H  =  u;  L  =  u; 
u=-2+u;  x=t.*(l-t);  v=l-2*t;  lamv=-u;  J=2;  H=-2+0*t; 
elseif  ell  >=  1/6 

half  1  =  f ind ( t <=  1 /2 )  ;  half  2  =  f ind ( t > 1/2 )  ; 

tl=t(halfl);  t2  =  l - 1 ( half  2 )  ; 

u=0*t ;  x=u;  v=u;  lamx=u;  lamv=u;  H=u; 

u(half 1 ) =-8* (l-3*ell) +24* (1 -4*  ell) *tl ; 

u(half 2) =-8* (1-3* ell) +24* (1 -4* ell) *t2 ; 

x (half  1 ) =  tl-4*(l-3*ell)*tl.  "2  +  4* ( 1 -4* ell ) *t 1  .  "3; 

x (half 2 ) =  t2-4*(l-3*ell)*t2 . " 2  +  4* ( 1 -4* ell ) *t2  .  "3; 

v (half  1 ) =  l-8*(l-3*ell)*tl  +  12*(l-4*ell)  .  *tl  .  "2; 

v (half 2 ) =-1  +  8* (1 -3* ell) *t2 -12* (1 -4* ell)  . *t2  .  "2 ; 

lamx (half  1 ) =24* ( 1 -4* ell )  ;  lamx (half  2 ) = -24* ( 1 -4* ell )  ; 
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lamv=-u;  J=2+6* ( 1 -4* ell ) ~ 2 ; 

Hnum  =  -8* ( 1 -6* ell ) “ 2 ;  H (half  1 ) =Hnum ;  H ( half  2 ) =Hnum ; 

else 

sl  =  find(t<  =  3*ell)  ;  s2=find(t>3*ell&t<(l-3*ell) )  ; 

s3=find(t>=(l-3*ell))  ;  s4  =  find(t>3*ell)  ; 

ta  =  l-  t(sl)/3/ell;  tb=(l-(l-t(s3)) /3/ell)  ; 

u ( s 1 ) = -2/3/ ell  * t a ;  u(s2)=0*s2;  u ( s3 ) = -2/3/ ell  * tb ; 

x ( si ) =ell * ( 1 -ta  .  ~3)  ;  x ( s2 ) =ell +0* s2  ;  x ( s3 ) = ell  * ( 1 - tb  .  ~3 )  ; 

v(sl)=ta.~2;  v(s2)=0*s2;  v ( s3 ) =-tb . " 2 ; 

lamx (sl)=2/9/ell.~2;  lamx (s4)=-2/9/ ell ~  2 ; 

lamv (si) =2/3/ ell  * ta ;  lamv (s4) =2/3/ ell*(l-(l-t (s4) ) /3/ ell)  ; 
J=4/9/ell;  Hnum=0;  H(sl)=Hnum;  H(s2)=Hnum;  H(s3)=Hnum; 

end 

results . cost  =  J  ;  results . solution . controls=u; 

results. solution. st at es=[x;v] ;  results. solution. H am iltoni an =H; 
results . solution . cost at es= [lamx; lamv] ;  results. solution. time=t; 

7, /  Plotting 

function  Plot _Re suit s ( result s ) 

7,  Assign  Constants 

J  =  results . cost  ;  soln  =  results . solution ; global  CONST , ell  =  C0NST . ell ; 
rBH  =  Bry  s  on_Ho  (  ell  )  ;  B  =  rBH.cost;  sBH  =  rBH  .  solut  ion  ;  7o  Bryson  and  Ho 
Np  =  length  ( soln)  ;  7»  soln  is  a  structure  array  containing  each  phase 
props={’defaultAxesFontSize  1  ,  ’ defaultLineLineWidth’  ,  .  .  . 

’  def  ault AxesXGrid  ’  ,  ’  def  ault AxesYGrid  ’  }  ;  ’/,  plot  properties 

old  =  get  (0  ,  props  )  ;  set  (0  ,  props  ,  { 12 , 2  ,  1  on  ’  ,  ’  on  ’  })  7»  reset  defaults 

%  Plot  States  and  Control 
figure(l) ,clf ,  mkr  =  { ’ bo  : ’ , ’ rd : ’ , ’ g~ : ’ } ; 
tb=sBH.time;  xb=sBH . states ;  ub=sBH . controls ; 

for  i=l:Np,  t=soln (i) . time ;  X=soln (i) . states ;  U=soln(i) . controls ; 
hl=subplot(3,l,l);  if  i==l,plot(tb,xb(l,:),,m-1),  end,  hold  on 
hi  =  plot (t , X ( :  ,  1 )  , mkr { i })  ;  ylabel(’x,)>  if  i ==1 ,  hleg  =  hl;  end 
h2=subplot(3,l,2);  if  i==l,plot(tb,xb(2,:),’m-,)I  end,  hold  on 
plot(t,X(: ,2) , mkr {i}),  ylabel(’v,)> 

h3=subplot(3,l,3);  if  i==l , plot ( tb , ub ( 1 , : ) , ’ m- ’ ) »  end,  hold  on 
plot(t,U(: ,1) , mkr {i}),  ylabel(’a,)J  hold  on 
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end,  xlabel  (  ’  time  1  )  7.  Compute  Cost  from  Bryson  &  Ho,  pg  .  122-123 

title  (hi  ,  sprintf  (  1  Cost_ {Bryson }  =  7, .  7  f  ,  Cost_{GP0CS}  =  7. .  7  f  ’  ,  B  ,  J  )  ) 
hbry  =  findobj  (hi ,  ’Color’  ,  ’magenta’)  ; 

hleg=legend ( [hbry , hleg] , ’  Bryson  &  Ho’,’  GPOCS ’,’ Location ’,’ S ’) ; 
hll  =  findobj  (hleg,  ’LineStyle ’  ,  ’  :  ’)  ;datl=get(hll , { ’ Xdat a  ’  ,  ’Ydata ’ })  ; 
hlm=findobj (hleg, ’Marker’ , ’o’) ;  dat2=get(hlm,{’Xdata’ , ’Ydata’}) ; 
hold(hleg,  ’on’)  ,  plot(hleg,datl-(l}(l)  ,datl{2}(l)  ,  ’bo’) 
plot (hleg, d at 2{1}(1)  ,dat2{2}(l)  ,  ’rd’  ,datl{l}(2)  ,datl{2}(2)  ,  ’ g~  ’  ) 
set(findall(gcf,  ’Col’,  [0,1,0]),  ’Color’,  [0,. 5,0])  7.  darker  green 
7.  Plot  Costates  and  Hamiltonian 

f igure (2)  , elf  ,  pb  =  sBH . costates  ;  Hb  =  sBH . Hamiltonian ; 
for  i=l:Np,  yl ={ ’ \ lambda_x ’ , ’ \ lambda_v ’ } ; 

t = soln ( i ) . t ime ;  L=soln(i) . costates ;  H= soln ( i ) . Hami It oni an ; 
hl=subplot(3,l,l);  if  i==l,plot(tb,pb(l,:),’m-’),  end,  hold  on 
hl=plot(t,L(:,l),mkr{i});  ylabel(yl{l}),if  i ==1 ,  hleg=hl ;  end 
h2=subplot(3,l,2);  if  i==l,plot(tb,pb(2,:),’m-’),  end,  hold  on 
plot(t,L(: ,2) , mkr{i}) , ylabel (yl{2}) 

h3=subplot(3,l,3);  if  i==l,plot(tb,Hb(l,:),’m-’),  end,  hold  on 
plot(t,H(:  ,1)  , mkr {  i  })  ,  ylabel(’H’)  ,  hold  on 
end,  xlabel(’time’),  title (hi, sprintf(’x_{max}  =  s’,rats(ell))) 
bry  =  findobj  (hi ,  ’Color  ’  ,  ’magenta’)  ; 

hleg=legend ( [bry , hleg] , ’  Bryson  &  Ho’,’  GPOCS ’,’ Location ’,’ NE ’) ; 
hll=findobj (hleg, ’LineStyle ’ , ’ : ’) ;datl=get(hll ,{’Xdata’ , ’Ydata’}) ; 
hlm=findobj (hleg, ’Marker’ , ’o’) ;  dat2=get(hlm,{’Xdata’ , ’Ydata’}) ; 
hold(hleg,  ’on’)  ,  plot(hleg,datl-(l}(l)  ,datl{2}(l)  ,  ’bo’) 
plot (hleg , dat 2 { 1 }  ( 1 )  , dat2{2} (1)  ,  ’ rd ’  , datl {1} (2)  , datl {2} (2)  ,  ’  g~  ’  ) 
set(findall(gcf,  ’Col ’,[0,1,0]),  ’Color’,  [0,. 5,0])  7o  darker  green 
7.  Plot  the  Path  Mutliplier  (if  a  path  constraint  is  assigned) 
if  isf ield (results . FUNCS ,’ path ’ )  &&  Np>=2 

f igure  (3)  ,plot(soln(2)  . time, soln  (2)  . multipliers ,paths(:,l),... 
’rd:  ’)  ,xlabel(’time’)  ,  ylabel(’ \mu ’),  title(’Path  Multiplier ’ ) 
end,  set  (0  ,  props  ,  old )  7.  return  to  previous  defaults 
7.7,  End  of  Breakwell  Problem 
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